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This  report  consists  of  two  code  manuals  and  an  article  recently  published  in  the 
proceedings  of  the  American  Nuclear  Society  Mathematics  and  Computation  Topical 
Meeting  held  in  Pittsburgh.  In  these  presentations,  deterministic  calculational 
methods  simulating  electron  transport  in  solids  are  detailed.  The  first  method 
presented  (Section  2)  is  for  the  solution  of  the  Spencer— Lewis  equation  in  which 
Electron  Motion  is  characterized  by  continuous  slowing  down  theory  and  a  pathlength 
formulation.  The  solution  to  the  standard  monoenergetic  transport  equation  foi 
electron  transport  with  isotropic  scattering  in  finite  media  is  given  in  Section  3. 
For  both  codes,  complete  flow  charts,  operational  instructions  and  sample  problems 
are  included.  Finally,  in  Section  4,  an  application  of  the  multigroup  formulation 
of  electron  transport  in  an  infinite  medium  is  used  to  verify  an  equivalent  S 
formulation.  For  this  case,  anisotropic  scattering  is  also  included. 
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1  Introduction 


The  final  report  to  follow  is  a  compilation  of  two  reports  and  an  article  recently  published 
in  the  proceedings  of  the  1991  American  Nuclear  Society  (A NS)  Topical  Meeting  of  the 
Mathematics  and  Computation  (M&C)  Division.  The  three  publications  deal  with  the 
development  of  deterministic  calculations!  methods  for  three  different  physical  characteri¬ 
zations  of  the  motion  of  electrons  in  a  solid  as  given  by  the  Boltzmann  transport  equation. 
These  methods  are  the  initial  step  in  developing  robust  algorithms  required  to  eventually 
predict  the  response  of  integrated  circuits  to  incident  electron  and  gamma  showers.  Highly 
accurate  numerical  solutions  have  been  generated  finding  use  as  benchmarks  to  which  al¬ 
gorithms,  developed  for  realistic  applications,  are  to  be  compared.  The  reports  not  only 
include  the  mathematical  theory  associated  with  the  solution  to  the  appropriate  transport 
equation,  but  also  the  computational  strategy  and  operational  instructions  for  the  accom¬ 
panying  computer  code.  A  computer  code  is  briefly  described  in  the  article  for  the  ANS 
Topical  Meeting  with  a  full  description  to  be  published  elsewhere. 

In  characterizing  the  motion  of  electrons  in  matter,  one  is  faced  with  a  spectrum  of  phys¬ 
ical  phenomena  all  requiring  a  model.  For  instance,  in  the  plasma  physics  application, 
collective  or  long  range  forces  play  a  dominate  role.  For  this  reason,  a  Vlasov  approxi¬ 
mation,  where  collisions  play  only  a  minor  role,  is  the  appropriate  form  of  the  transport 
equation.  In  the  solid  state  regime,  where  collisions  provide  the  forces  for  electron  motion 
at  high  energy  and  the  forces  from  fields  are  relatively  weak,  a  linear  transport  equation 
provides  an  adequate  description.  In  this  formulation,  the  major  assumption  is  that  the 
atomic  electric  fields  are  localized  (within  the  Debye  sphere)  allowing  the  application  of 
the  collision/streaming  statistical  model  given  by  the  linear  Boltzmann  equation.  There 
exist  further  approximations  which  lead  to  several  simplified,  but  useful,  models  such  as: 


•  the  Spencer- Lewis  (S-L)  equations, 

•  the  monoenergetic  transport  equation,  and 

•  the  multigroup  transport  equation. 


These  are  the  three  specific  areas  detailed  in  the  following  chapters  in  this  report. 

In  the  S-L  formulation  (in  §2),  the  continuous  slowing  down  approximation  (CSDA)  is 
enforced  where  electrons  are  restricted  to  a  specific  energy  loss  for  a  given  penetration  dis¬ 
tance  or  pathlength.  When  the  interaction  probabilities  are  assumed  constant,  in  addition, 
the  solution  to  the  S-L  equation  can  be  obtained  from  a  combination  of  a  multiple  collision 
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expansion  and  a  moments  reconstruction.  The  analytical  representation  in  terms  of  path- 
length  can  be  evaluated  to  an  exceptionally  high  accuracy  to  establish  a  true  benchmark. 

The  monoenergetic  transport  equation  in  a  finite  medium  (in  §3)  is  solved  with  the  highly 
accurate  Fjv  method.  The  algorithm  features  heterogeneous  media  with  a  uniform  isotrop¬ 
ically  emitting  source  in  each  slab.  Currently,  this  algorithm  is  limited  to  isotropic  scat¬ 
tering  with  the  intent  to  include  anisotropic  scattering  and  down  scattering  in  the  more 
advance  multigroup  version. 

The  final  presentation  (in  §4)  concerns  the  multigroup  formulation  of  electron  down  scat¬ 
tering  in  an  infinite  medium.  The  solution  with  anisotropic  scattering  is  obtained  via  a 
numerical  Fourier  transform.  This  accurate  multigroup  benchmark  is  compared  to  a  stan¬ 
dard  discrete  ordinates  approximation  from  which  the  order  of  the  error  for  various  spatial 
approximations  relating  the  average  to  the  cell-edge  flux  is  confirmed. 
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SLEET 


The  ^pencer- Lewis  Equation  of  Electron  Transport:  A  Code  Manual 


B.  Ganapol* 

Department  of  Nuclear  and 
Energy  Engineering 

University  of  Arizona 
Tucson,  Arizona 


Abstract 

The  Spencer-Lewis  equation  of  electron  transport  with  energy  independent  cross  sections  and  a  screened- 
Rutherford  scattering  kernel  in  an  infinite  medium  is  solved  using  moments  reconstruction  coupled  with  the 
multiple  collision  formulation.  The  solution,  evaluated  to  a  high  degree  of  accuracy,  finds  use  as  a 
benchmark  against  which  purely  numerically  based  algorithms  can  be  tested.  The  evaluation  features 
general  anisotropic  scattering  for  a  localized  source  or  for  a  uniformly  distributed  source  in  a  half-space. 
The  operation  of  the  SLEET  code  which  embodies  the  numerical  solution  is  described  in  detail  in  this 
document. 


t  Work  performed  for  Rome  Air  Development  Center  at  Hanscom  AFB 
CALSPAN* 

Contract  No.  F30602-88-D-0026,  S-0-7542 
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I.  INTRODUCTION 


With  the  ever  decreasing  dimensions  of  integrated  circuits  and  the  necessity  for  shielding  these 
sensitive  electronic  components  against  penetrating  electrons,  a  renewed  interest  in  electron  transport  theory 
has  recently  been  inspired.  Because  of  the  complex  nature  of  electron  transport  in  solids,  direct  numerical 
simulation  of  electron  motion  by  Monte  Carlo  calculations  has  found  widespread  use.  Monte  Carlo 
simulations  are  flexible  and  can  incorporate  a  variety  of  physical  processes  in  complicated  three-dimensional 
configurations.  At  the  same  time  however,  probabilistic  analyses  can  require  excessive  computational  effort 
and  therefore,  may  not  be  the  most  efficient  calculational  tool  for  all  problems.  Thus,  new  approaches 
leading  to  a  numerical  characterization  of  electron  transport  in  solids  have  and  will  continue  to  be 
developed. 

In  a  classic  paper,  Lewis1  proposed  a  Boltzmann  equation  description  of  electron  motion  in  terms  of 
path  length.  In  this  formulation,  angular  deflection  is  specified  by  a  scattering  term  representing  deflection 
without  energy  loss;  and  energy  degradation  is  characterized  by  a  prescribed  range-energy  relation.  Lewis's 
description  assumed  that  the  electrons  continuously  lose  energy  with  motion  thus  the  statistical  nature  of  the 
slowing  down  process  is  neglected.  The  resulting  Boltzmann  equation,  called  the  Spencer-Lewis  equation, 
was  first  solved  numerically  by  Spencer2  for  an  infinite  homogeneous  medium  by  using  moments  to 
reconstruct  the  electron  dose  from  a  functional  fit  based  on  physical  augments.  In  addition  to  the  moment 
function  fitting  procedure  of  Spencer,  a  significant  effort  has  been  given  to  the  application  of  SN  methods 
to  these  problems.2 

The  development  of  SN  and  other  numerical  methods  for  use  in  electron  transport  applications  has 
given  rise  for  the  need  to  assess  the  accuracy  of  proposed  numerical  algorithms  and  to  provide  the  assurance 
that  these  algorithms  have  been  properly  programmed.  Analytical  solutions  for  representative  electron 
transport  problems  which  can  be  evaluated  to  a  high  degree  of  accuracy  serve  as  standards  or  benchmarks 
for  comparison  with  the  approximate  solutions.  Such  a  solution  has  been  developed  for  the  S-L  equation 
characterizing  electron  transport  in  a  simplified  physical  setting.  Specifically,  electrons  are  assumed  to  slow 
down  in  an  infinite  medium  after  emission.  Furthermore,  all  scattering  and  absorption  interactions  are 
assumed  to  be  energy  independent  (except  when  otherwise  stated)  with  the  kinematics  of  the  scattering 
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process  described  by  a  screened-Rutherford  law. 


The  method  of  numerical  solution  employs  a  flux  reconstruction  technique  a  la  Spencer,  however, 
the  functional  form  of  the  desired  flux  distribution  is  not  known  a  priori.  The  reconstruction  is  based  on 
an  expansion  in  Legendre  polynomials  coupled  with  the  multiple  collision  formulation.  In  this  way,  the 
solution  can  be  evaluated  to  a  high  degree  of  accuracy  (4-5  digits).  The  SLEET  (Spencer-Lewis  Equation 
of  Electron  Transport)  program  has  been  written  specifically  for  this  purpose.  This  report  is  a  manual 
describing  the  SLEET  code  in  detail  including  theory  and  operation. 


II.  THE  TRANSPORT  DESCRIPTION 
A.  The  Spencer-Lewis  (S-L)  Equation 

The  S-L  equation  describing  the  motion  of  electrons  experiencing  deflection  and  energy  loss  in  one¬ 


dimensional  plane  homogeneous  medium  with  azimuthal  symmetry  is 

r  i  f2’  r1 

i  * "  Jr  *  f<x-',-!|  •  W)  J0  M'  j.™  p('‘«-s)  f(XJ,'-s)  * Q(w) 


where 


s  a  cumulative  electron  path  length 


x  a  distance  measured  along  coordinate  axis 
f  a  angular  density  of  electrons 
X  a  total  mean  free  path 

p(/i0,s)d/x0  *  probability  of  an  elastic  scattering  event  producing  a  deflection  between  n0  and  n0  +  d/r0 
<t>  a  azimuthal  angle 

q  dx  d/i  ds  •  volumetric  rate  of  external  injection  of  electrons  at  position  x  in  direction  n  with  path 
length  s. 

In  the  customary  fashion,  the  scattering  kernel  is  assumed  to  be  expandable  in  terms  of  Legendre 

polynomials  P(. 


pOvs> 


(2) 


oo 

Z  ^(s)  p^<>>  •  wo  - 1 


and  Eq.  (1)  becomes 


k +  **  It  +  A<b  f(x’M’s)  =  j  8(^’s)  f(x^s)  +  Q(w) 


where 


oo 

/— ^,S)  s  ^  W£(S)  P till)  P i(n')  . 


This  equation  will  be  solved  in  an  infinite  medium  with  constant  properties  and  with  the  following 

boundary  condition  (for  a  localized  source); 

lim  f(x,/i,s)  ■  0 
|x|— *oo  . 

In  addition,  electrons  are  physically  restricted  to  positive  path  lengths  requiring 
f(x,/i,s)  *  0,  s  <  0  . 


B.  Extended  Transport  Correction 

Typically  in  electron  transport  applications,  the  series  representation  for  the  scattering  kernel  given 
by  Eq.  (3b)  requires  many  terms  to  provide  an  adequate  description  of  the  scattering  event.  An  extended 
transport  correction4  has  been  proposed  as  a  way  to  reduce  the  number  of  terms  in  the  scattering 
approximation  while  maintaining  the  desired  accuracy.  In  particular,  the  series  for  is  replaced  by  a 

series  with  modified  coefficients  ui\  truncated  at  L  plus  a  delta  function  to  represent  the  forward  directed 
component;  thus  Eq.  (3a)  becomes 

Js  +  *  dx  *  pj  *  jr  |  d^'  g*0*'-V»)  +  q(x^,s)  (4a) 

with  the  extended  transport  corrected  kernel  and  mean  free  path  given  by 


sV-m)  a  Z  ^T~  W) 


u)t  «  K-<^+1]/[l-^+i]  ,  0  <  l  <  L 
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(4d) 


C.  Screened-Rutherford  Scattering  Kernel 

Since  electrons  are  negatively  charged  particles,  the>  are  deflected  by  the  electric  fields  associated 
with  the  electronic  cloud  as  well  as  the  nucleus  of  the  host  medium  atoms.  For  this  situation,  the  screened- 
Rutherford  scattering  kernel  is  an  appropriate  approximation  with  the  following  Legendre  coefficients 
derived  in  Appendix  A: 

w0  “  I  (5a) 

w,  =  1  +  r?-r?(  1+77/2)  £n(l+2/*7)  (5b) 

=  (2£+l)  (I+77)  u^-(£+l)u^_i  (5c) 


and 

n  =  3.2  x  IO-3  z2'3  E(keV)  (5d) 

1  =  Z(Z+|)  .  _J_ 

A  A  4  £(MeV)2  1+^/2 

where 


Ee  *  electron  rest  mass  energy  (0.511  MeV) 
E  =  average  electron  energy 
Na  =  Avogodro's  number 
A  =  atomic  weight  of  the  scattering  center 
Z  =  atomic  number  of  scattering  center 
r,  «  electron  radius  (2.8179  x  10'13  cm). 


D.  Source  Distribution 

The  source  is  assumed  to  be  located  at  the  zero  plane  (x=0)  and  to  emit  electrons  having  a  zero  path 
length  with  an  angular  distribution  specified  by  Q (71),  thus 
q(x,A*,s)  =  Q(m)  «(x)  6(s). 

The  angular  source  distribution  will  be  of  two  types 


Ab 


1/2,  isotropic  emission 


\^Q(u-p0),  azimuthally  symmetric  beam  emission 
where  n0  is  the  cosine  of  the  angle  of  incidence  of  the  source  electrons. 


HI.  RECONSTRUCTION  OF  THE  SCALAR  DENSITY 
A.  Multiple  Collision  Formulation 

By  decomposition  of  f(x,/x,s)  into  its  collisional  components  fn,  Eqs.  (4a),  (3c)  and  (3d)  are 
transformed  into  the  following  multiple  collision  representation: 


IL  f0(x,|i,s)  =  Q (m)  6(x)  6(s) 


IL  fn(x,p.s) 


! 

dp'  g’(/r'- 

J-i 


•p)  fn- i  (x,p  ,s)  ,  n  =  1,  2,  .  .  . 


lim  fn(x,p,s)  =  0 
I  x  l—oo 


with  the  total  angular  flux  given  by 


OO 

(x,/i,s)  =  ^  fn(x,M,s)  . 


By  application  of  integral  transport  theory  and  the  definition  of  the  similarity  variable  rj  =  x/s,  the  scalar 
flux  solutions  to  Eqs.  (7)  for  n  *  0,1  are 


f0(x,s)  *  -  0O(» ,)  $d-|n|) 

where 


t&ofa)  5  Q<7)  . 


and 


f\(x,s)  =  e-*/A  01(»7)  0(1 -hi) 


(8c) 


with 


V'xto)  S  f  d/*  f  d r\  g'(ri'->n)  +  j  d(i  f  dr?'  g *(r?'-.?i)  • 

J  - 1  Jr?  « i?  J  - 1 


(8d) 


0(u)  is  the  Heaviside  step  function.  From  the  correspondence  with  the  time-dependent  neutron  transport 
case5,  it  can  be  shown  that 

fn(x,M,s)=  Fn(M,r?)<Kl-H)  (9) 

where  Fn  satisfies  the  reduced  collision  equation 

,1 


r„(M) 


I 

n  d/i' 
•  -1 


n)  Fn.jC/i'.r?) 


for  n  >  1  with  the  condition 


lim  Fn(/i,r?)  =  0 

hl-i 

Thus  the  scalar  density  to  be  evaluated  is  given  by 
f(x,s)=  ^n(r?)  0(1  -|r?|) 


(10a) 


/10b) 


(Ha) 


n=0 


with 


=  |  d/i'  Fn(/i',r?)  . 


(lib) 


B.  Legendre  Series  Expansion  for  t/)n(r?) 

Since  the  collisional  process,  even  for  highly  forward  peaked  scattering,  provides  a  mechanism  for 
redistribution  of  particles  in  phase  space  and  therefore  a  smoothing  of  highly  directed  sources,  Fn  as  a 
function  of  /j  and  r?  (for  n  sufficiently  large)  is  expected  to  be  continuous  and  uniformly  bounded.  By 
taking  this  behavior  into  account  and  since  Fn  is  nonzero  for  |r?|  <  1,  the  following  expansion  of  Fn  in 
Legendre  polynomials  in  the  variable  r?  seems  entirely  appropriate: 


Ao 


F„  (/!,>?) 


(12a) 


OO 

^  ^ifn  k(#i)  Pk 


<t?) 


k=0 


where 


fn,k(M)  s  |  dr?  Pk(r?)  F n(n,tj)  ; 

and  for  the  scalar  component 
oo 

Mv)  =  Y,  ^T~  f°n'k  W  • 

k=0 

where  the  coefficients  f°n  k  are 

fon,k  =  |  dM  fn.k(»7)  *n0l>  • 

The  expansion  coefficients  are  found  directly  from  Eq.  (10a)  by  substitution  of  Eq.  (12a), 

Pk(r?)  P i(n)  and  integrating  over  r?  [-1,  1]  and  /j  [-1,  1]  to  yield  the  recursion  relation 

,1 


C  *  W  | 


dr?  Q(r?)  P^(r?) 


fo,k  «  j  d»»  P lW  Pk(r?)  Q(r?) 
and  for  k  >  1,  n  >  1 

(k+n)  «  -  (k-n-I)  f*,k.2  ♦  |^f  [t  fft.j  ♦  (M)  fJjILj  + 

+  n  w/[fn-l,k  “  fn-l.k-2] 


with 


*  0  . 


C.  Spatially  Distributed  Source 

The  scalar  density  h(x,s)  for  a  spatially  distributed  source  of  the  form 
q(x,/i,s)  -  U(x)  Q(m)  6( s) 


(12b) 

(13a) 

(13b) 

multiplying  by 

(14a) 

(14b) 

(14c) 

(14d) 

(15a) 
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can  be  represented  by  the  convolution  integral 


r 

h(x,$)  =  c 

J  -OC 


dx'  U(x-x')  f(x',s) 


(15b) 


where  f(x,s)  is  given  by  Eq.  (11a).  Thus  Eq.  (15b)  becomes,  after  substitution  of  Eqs.  (11a)  and  (13a), 


oo  oo 

h(x,s)  e‘*A*  X  (S/A*)n  f"-k  IkM 


(16a) 


k-0 


n=0 


where 


Ik(t|)  *  J  dtf 


U(s-xn')  Pk(i)')  . 


(16b) 


D.  Power  Law  Variation  of  A*'1 

If  A*-1  is  assumed  to  be  of  the  form 

A**i(s)  -  A0(s/s0)*  .  -1  <  P  <  oo  ,  (17) 


then  the  solutions  obtained  thus  far  can  be  shown  to  have  the  same  representation  with  the  following 
replacement 


and 


(18a) 


[k+n(l+0)]  -  -  [k-n(l+0)-IJ  {*^{0)  ♦  0?)  ♦  (*♦!)  fft-i  (fll  ♦ 

♦  n  «£  f<tlik  (*>  -  fn-lik-2  W]  (18b) 

fo,k  (£)'*  fo,k  .  (18c) 
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E.  Uncollided  and  First  Collided  Fluxes 


As  a  result  of  the  localized  nature  of  the  source  in  x  and  s,  the  uncollided  scalar  flux  and  possibly 
the  first  collided  flux  are  singular.  Since  the  Legendre  expansion  converges  slowly  near  singularities,  it  is 
advantageous  to  remove  these  components  from  the  expansion.  This  is  done  by  obtaining  the  uncollided 
and  first  collided  fluxes  analytically  from  the  solution  of  Eqs.  (7a)  and  (7b)  for  n  =  1  (for  a  beam  source)  to 


give  for  an  isotropic  source  [Q(/i)  =  Q0/2] 

f0(x,s)-  <'*» 

and  for  a  beam  source 

f0(x,s)  =  Q0  Si  T)-n0)  (19b) 

f  ,(x,s)  =  Q0  [8(n0-n  )  J-(r?)  +  fKn-n 0)  J+(rj)J  (19c) 

L 

JHv)  s  ^  ^  P*(Mo)  Vf  09d) 

t=0 

Yq  =  ln(l±/i)/(/i0-»?)|  (I9e) 

Y;  =  ±(Ur?)  +  n0  Y0*  090 

Y |  -  /t0  Yf.,  -  ^  Y U  -  \  [P*(t?)  -  P*.2(t?)]  .  (19g) 

For  a  uniformly  distributed  isotropic  source  for  x  <  0 
%-»/**  ,  17  <  -  1 

f0(x,s)  =  Q0  ■  e*,/r  (*-»l)/2  ,  H  S  1  (20a) 

0  ,  v  >  I 


and  for  a  beam  source 

r ,  ?7  <  —  1 

f0(x,s)=  Q0  ■  e-»/A*  <7(ji0-r?)  ,  hi  <  1  (20b) 

0  ,  rj  >  1 
v. 

IV.  Numerical  Evaluation 


The  expression  for  the  scalar  density  obtained  for  a  localized  source  at  x  ■  0  is 


n, - 1  oo 

C,s)  =  2]  fn<X’S>  +  X  ^T1  f“(S)  W 


where  the  discontinuous  uncollided  (n  =  0)  and  first  collided  (n  =  1)  densities  have  been  removed  from  the 
Legendre  expansion  in  order  to  expedite  convergence  and 


A  further  simplification  can  be  realized  by  noti'*'’  that  for  n,  =  2 


f0(s)  =  i  M0(s) 


f°  *  A.  f0 
"  *0,0  +  )*  *i,( 


where 


M0(s)  *  dx  f(x,s)  . 


From  particle  conservation  M0(s)  ■  1  ,  and  from  Eq.  (14a) 


fS.o  »  1  . 


f0(s)  -  1  -  (1  +  jf)  /s 


for  a  beam  source  and 

f0(s)  *  (1  -  e**/A‘)/s 

for  an  isotropic  source.  The  final  expression  to  be  evaluated  is  therefore 


n,  - 1  00 

t,s)  -  2^  fn(s)  ♦  fk(s)  Pk(n)  • 


To  evaluate  this  expression  the  following  numerical  operations  are  required: 
recursive  evaluation  of  f|J  k 
multiple  collision  series  summation  (sum  over  n) 
moments  reconstruction  (sum  over  k)  . 
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A.  Expansion  Coefficients 


The  major  numerical  error  associated  with  the  recursion  relation  for  f£k  is  round-off.  To 
determine  f£k,  the  coefficients  for 
n  =  0,  1,  .  .  .  Ln 

k  =  0,  1 - Lk 

t  =  0,  1,  . . .  lk 

are  required  due  to  the  telescoping  nature  of  the  recurrence  relation.  LN  and  LK  are  the  maximum  number 
of  terms  allowed  in  the  multiple  collision  and  Legendre  series  evaluation  respectively. 

B.  Convergence  Iteration 

The  infinite  series  are  evaluated  using  a  simple  engineering  estimate  where  the  series  is  considered 
converged  when  each  of  IEK  terms  terms  produces  a  relative  error  of  less  than  a  specified  amount  (ERN 
and  ERK  for  the  multiple  collision  and  Legendre  series  respectively).  The  multiple  collision  series  is 
evaluated  first  (for  each  x  and  s)  for  a  truncation  error  ERN  »  ERK./10.  If  the  series  converges,  a  further 
test  is  performed  to  determine  whether  the  multiple  collision  summation  (over  n)  has  been  evaluated  with 
sufficient  digits  to  give  the  desired  accuracy.  This  is  accomplished  by  comparison  of  each  term  of  the 
Legendre  sum  with  the  final  result.  If  the  final  outcome  has  been  obtained  as  a  result  of  addition  and 
subtraction  of  much  larger  numbers  than  the  final  result,  the  flux  may  largely  be  noise  if  the  individual 
terms  did  not  contain  enough  digits.  If  true,  a  second  iteration  begins  where  the  multiple  collision  sum  is 
recalculated  with  a  reduced  truncation  error  ERN/10r  ,  where  r  is  the  estimated  number  of  digits  required, 
and  the  Legendre  series  is  again  summed.  The  process  is  continued  until  ERN  <  10'20  or  the  flux  changes 
less  than  ERK.  between  successive  iterations. 

C.  Convergence  Acceleration 

For  highly  anisotropic  scattering,  the  Legendre  series  can  be  slowly  convergent  or  even  numerically 
divergent.  The  divergence  results  from  round-off  error  that  produces  inaccuracies  in  the  terms  of  the 
summation  eventually  leading  to  numerical  noise  and  subsequent  divergence.  For  isotropic  and  beam 
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Fig.  1  SLEET  logic  flow  chart 


sources  (at  x*0),  the  convergence  of  the  Legendre  series  can  be  accelerated  by  removal  of  the  singularities 

at  tj  -  ±  1.  This  is  accomplished  by  representation  of  the  flux  that  has  had  at  least  are  collision  by 

oo 

(I -ij*)"»  [f(x,s)  -  f0(x,s)J  =  A^(k)  Pk(r?)  (25a) 

k=0 

The  motivation  for  this  expansion  is  that  the  nth  derivative  of  rl>^  possesses  singularities  at  r?  =  +  1  and  the 
above  expansion  tends  to  remove  these  singularities.  The  expansion  coefficients  are  obtained  from  the 
recursion  relation: 

a[[0)(s)  =  fk(s) 

A^m)(s)  =  ( 1  -Tf2k  >  Ak"1'11  '  1*  A^j11  -  73k  A&*1*  (25b) 

llk  h  (k+lKk+2)/ [2k+l)(2K+3) 

— iJsiL)! -  +  _ k! _ 

l2k  (2k+lX2k+3)  (2k-l)(2k+l) 

73k  =  (k+l)k//(2k-lX2k+l)  . 

Because  the  function  expanded  is  smoother  than  f(x,s)  itself,  the  series  in  Eq.  (22)  will,  in  general, 
converge  more  rapidly  than  the  series  in  Eq.  (21a). 

V.  The  SLEET  Code 

The  algorithm  described  above  was  coded  in  a  FORTRAN  77  program  called  SLEET  which  is  an 
acronym  for  the  jipencer-Lewis  Equation  of  Electron  Transport.  A  flow  diagram  indicating  the  logic  of  the 
program  is  given  in  Fig.  1.  Each  cell  has  been  given  a  number  indicating  the  calling  sequence  allowing  the 
original  call  to  be  traced.  The  program  is  composed  of  IS  routines  each  with  a  specific  function 
summarized  in  the  flow  chart. 

A.  Input  Options 

The  SLEET  input  description  is  given  in  Fig.  2.  the  input  file  is  called  SLEET.DAT  and  must  exist 
at  execution.  Free  format  is  used  where  the  entries  on  each  line  must  be  separated  by  at  least  a  space. 

The  following  options  are  available: 

1 )  Isotropic  or  beam  source  angular  distributions. 

2)  Delta  function  or  uniform  half-space  (x<0)  spatial  source  distributions 
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Fig.  2 


SLEET  input 


**************************************************************** 
INPUT  DESCRIPTION 

CARD  1  (18A4) 

IDEN  -  PROBLEM  IDENTIFICATION 
CARD  2 

IS  -  1  ISOTROPIC  SOURCE 

-  2  ANISOTROPIC  SOURCE 

JTI  -  0  DELTA  FUNCTION  SOURCE  AT 
X-0 

-  1  CONSTANT  SOURCE  IN  LEFT 

HALF-SPACE ( X< 0 ) 

QO  -  SOURCE  STRENGTH 

AMO  -  COSINE  OF  ANGLE  OF  INCIDENCE  OF  SOURCE 
ELECTRONS 

CARD  3 

LL  -  ORDER  OF  SCATTERING  KERNEL  APPROXIMATION 

ITC  —1  SCREENED-RUTHERFORD  SCATTERING  KERNEL 

-  1  READ  IN  SCATTERING  COEFFICIENTS 

ITRN  — 1  USE  TRANSPORT  CORRECTED  CROSS  SECTIONS 

-  1  DO  NOT  USE  TRANSPORT  CORRECTED  CROSS 

SECTIONS 

CARD  4 A  IF  ITC-1 
ALA  -  TOTAL  MEAN  FREE  PATH 

CARD  4B  IF  ITC— 1 

ETA  -  SCREENING  PARAMETER  FOR  SCREENED 

RUTHERFORD  SCATTERING  CROSS  SECTION 

ALA  -  TOTAL  MEAN  FREE  PATH 


CARD  5 

BE  -  POWER  OF  PATH  DEPENDENCE  OF  NUMBER 
OF  SECONDARIES  ( ALA( S ) - ( S/SlO ) **BE ) 
SIO  -  CROSS  SECTION  PARAMETER 

CARD  6 

L S  -  NUMBER  OF  PATH  GRID  POINTS 
LX  -  NUMBER  OF  SPACE  GRID  POINTS 

JT  —1  READ  IN  PATH  GRID  POINTS 

-  1  READ  IN  PATH  MESH  SPACING 

JX  —1  READ  IN  SPACE  GRID  POINTS 

-  1  READ  IN  SPACE  MESH  SPACING 
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CARD  7 

IDFF— 1  COMPARISON  WITH  DIFFUSION  THEORY 
(BE-0  ONLY) 

-  1  SCALAR  FLUX  CALCULATION 
CARD  8 

LN  -  MAXIMUM  NUMBER  OF  TERMS  IN 

MULTIPLE  COLLISION  EXPANSION 
OF  MOMENTS 

LK  -  MAXIMUM  NUMBER  OF  TERMS  IN 
LEGENDRE  SERIES  EXPANSION 

IEK  -  NUMBER  OF  CONSECUTIVE  TERMS 
SATISFYING  RELATIVE  ERROR 
CRITERION  IN  LEGENDRE  EXP- 
SION( 3) 

NOTE : I F  IEK.LT.0  SUPPRESS  ANALYTIC  DETERMINATION  OF 
OF  FIRST  COLLIDED  FLUX  WHEN  APPROPIATE 

LM  -  NUMBER  OF  CYCLES  FOR  ACCELERATED 
CONVERGENCE ( 3 ) 

CARD  9 

ERK  -  RELATIVE  TRUNCATION  ERROR  FOR 
LEGENDRE  EXPANSION  1 . 0E-04 ) 

CARD  10 

ITR  — 1  PRINT  OUTPUT  ON  PRINTER  (132  COLUMNS) 

-  1  PRINT  OUTPUT  ON  TERMINAL  (72  COLUMNS) 

IP  —1  DO  NOT  PRINT  INPUT  EDIT 

-  0  PRINT  INPUT  EDIT  WITH  OUTPUT  EDIT( TAPE21 ) 

-  1  PRINT  INPUT  EDIT  ON  TAPE 2 3 

IXY  —1  S  ON  ABCISSA, X  ON  ORDINATE  IN  EDIT 
-1  X  ON  ABCISSA, S  ON  ORDINATE  IN  EDIT 

I  PL  —1  NO  PLOTS  REQUESTED 

-  1  WRITE  PLOT  FILE(TAPE22) 

CARD  11 

ITPE  —1  READ  FNK  FROM  TAPE24  AND  COMPLETE 
CALCULATION 

-  0  CREATE  FNK  ON  TAPE24  AND  COMPLETE 

CALCULATION 

-  1  CREATE  FNK  ON  TAPE24  AND  STOP 
CARD  12  IF  ITC-1 

FTC  -  SCATTERING  COEFFICIENTS  (IC-0,LL) 

CARD  13A  IF  JT— 1 
SS  (  IT ) -  PATH  GRID  POINTS ( IT-1 , LS ) 

CARD  13B  IF  JT-1 

SO  -  INITIAL  PATH  LENGTH 
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DS  -  PATH  INTERVAL 
CARD  14A  IF  JX— 1 

XX(IX)-  SPACE  GRID  POINTS  (IX-1,LX) 

CARD  14B  IF  JX-1 

XO  -  INITIAL  SPACE  GRID  POINT 
DX  -  SPACE  MESH  SPACING 
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3)  Screened-Rutherford  transport  corrected  scattering  coefficients 

4)  Power  law  variation  of  A*-1 

5)  Comparison  with  diffusion  theory. 
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1)  Problem  1 


The  first  problem  is  the  basic  isotropic  scattering  and  source  problem  where  the  source  is  at  x  •  0. 
The  flux  is  determined  at  x  -  1(1)5  for  paths  S  -  1(1)20. 

Along  with  the  fluxes,  diagnostic  information  indicating  the  convergence  properties  of  the 
reconstruction  at  each  time  and  space  point  is  given  in  the  output  edit.  The  information  is  contained  in 
coded  number  under  the  title  "NUMBER  OF  TERMS  FOR  CONVERGENCE".  Each  entry  corresponds  to 
the  x  and  s  in  the  table  above.  If  no  infinite  series  evaluations  were  performed  for  a  particular  x,t  point, 
then  a  zero  entry  is  found.  This  will  be  the  case  in  the  forbidden  region  (s<x)  or  if  x  and  s  are  equal  and 
the  flux  is  the  uncollided  contribution  only.  Otherwise  a  number  of  the  form  bxx.yyyuEizz 
will  be  found.  The  number  is  coded  as  follows: 
u  is  the  number  of  convergence  accelerations  required 

yyy  is  the  maximum  number  of  terms  required  for  convergence  of  the  multiple  collision  series 

-zz  indicates  the  Legendre  series  has  converged  to  the  desired  relative  error(ERK./10)  with  zz 

being  the  estimated  relative  error  exponent  for  the  multiple  collision  series 
+zz  indicates  the  Legendre  series  did  not  converge  to  the  desired  relative  error  with  zz  now 

representing  the  exponent  of  the  minimum  relative  error  for  the  Legendre  series 
b«-  indicates  that  more  than  one  convergence  iteration  (IV.B)  has  been  performed  in  order  to 

achieve  the  desired  precision 

b-blank  indicates  that  no  convergence  iterations  were  necessary  to  achieve  the  desired  precision 

xx  represents  the  number  of  terms  required  for  convergence  of  the  Legendre  series  unless  a  + 

sign  precedes  zz.  In  this  case,  the  series  did  not  converge  and  xx  is  the  term  in  the  series 
producing  the  least  relative  error.  The  flux  has  been  set  equal  to  the  partial  sum  up  to  and 
including  this  term. 
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a)  Input:  File  SLEET . DAT 

SAMPLE  PROBLEM  1: ISOTROPIC  SOURCE  AT  X-0 , ISOTROPIC  MEDIUM 
1  0  1.0  1.0 
Oil 
1.0 

0.0  1.0 
20  5  1  1 
1 

90  60  3  3 
1 .0E-04 
111-1 
0 

1.0 

1.0  1.0 
1.0  1.0 
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b)  Input  Echo:  File  03 


SAMPLE  PROBLEM  1: ISOTROPIC  SOURCE  AT  X-0 , I SOTROPIC  MEDIUM 


IS  - 

1 

1- ISOTROPIC  SOURCE 

2 - AN I SOTROPIC  SOURCE 

JTI- 

0 

0-DELTA  FUNCTION  SOURCE  AT  X-0 

1-CONSTANT  SOURCE  IN  LEFT  HALF-SPACE ( X<0 

Q0  - 

1.0000E+00 

SOURCE  STRENGTH 

ALA  - 

1.0000E+00 

MEAN  FREE  PATH 

LL  - 

0 

ORDER  OF  SCATTERING  APPROXIMATION 

IT C  - 

1 

-1-USE  TEST  KERNEL 

1-READ  COEFFICIENTS 

ITRN- 

1 

-1-USE  TRANSPORT  CORRECTED  CROSS  SECTIONS 
1-DO  NOT  USE  TRANSPORT  CORRECTION 

SCATTERING  COEFFICIENTS( L-0-  0) 
1.0000E+00 


BE 

* 

0.0000E+00 

EXPONENT  IN  VARIATION  OF 
( C ( S ) -C* ( S/SI  0 ) *  *BE ) 

SI0 

m 

1 . 0000E+00 

PARAMETER  IN  VARIABLE  C 

LS 

m 

20 

NUMBER  OF  PATH  POINTS 

LX 

m 

5 

NUMBER  OF  SPACE  POINTS 

JT 

- 

1 

-1-READ  IN  PATH  POINTS 
1-READ  IN  MESH  SPACING 

JX 

* 

1 

-1-READ  IN  SPACE  POINTS 
1-READ  IN  MESH  SPACING 

SO 

m 

1 . 0000E+00 

INITIAL  PATH 

DS 

m 

1 .0000E+00 

PATH  INTERVAL 

XO 

m 

1 . 0000E+00 

FIRST  POSITION 

DX 

1 . 0000E+00 

SPACE  INTERVAL 

PATH 

POINTS 

1 . 0000E+00  2 . 0000E+00  3.0000E+00  4.0000E+00  5.0000E+00 

6. 0000 E+ 00  7 . 0000E+00  8.0000E+00  9.0000E+00  1.0000E+01 

1 . 1000E+01  1 . 2000E+01  1.3000E+01  1.4000E+01  1.5000E+01 

1 . 6000E+01  1 . 7000E+01  1.8000E+01  1.9000E+01  2.0000E+01 

SPACE  POINTS 

1 . 0000E+00  2 . 0000E+00  3.0000E+00  4.0000E+00  5.0000E+00 

IDFF-  1  -1-DIFFUSION  COMPARISON 

1-SCALAR  FLUX  CALCULATION 

LN  -  90  MAXIMUM  NUMBER  OF  TERMS  IN 

MULTIPLE  COLLISION  EXPANSION 
ERN  -  1 . 0000E-05  RELATIVE  TRUNCATION  ERROR 

FOR  MULTIPLE  COLLISION  EXPANSION 
IEN  -  3  NUMBER  OF  TERMS  SATISFYING  RELATIVE 


A22 


ERROR  CRITERION 


LK  -  60  MAXIMUM  NUMBER  OF  TERMS  IN 

LEGENDRE  EXPANSION 

ERK  -  1.0000E-04  RELATIVE  TRUNCATION  ERROR 

IEK  -  3  NUMBER  OF  TERMS  SATISFYING  RELATIVE 

ERROR  CRITERION 

LM  -  3  NUMBER  OF  PA  ES  FOR  CONVERGENCE ( JTI-0 ) 

ITR  -  1  -1-OUTPUT  TO  PRINTER(132  COLUMNS) 

1-OUTPUT  TO  TERMINAL( 72  COLUMNS) 

IP  -  1  -1-DO  NOT  PRINT  INPUT  EDIT 

0-PRINT  INPUT  EDIT  ON  TAPE21 
1-PRINT  INPUT  EDIT  ON  TAPE23 

IXY  -  1  -1-S  ON  ABCISSA  AND  X  ON 

ORDINATE  IN  OUTPUT  EDIT 
1-X  ON  ABCISSA  AND  S  ON 
ORDINATE  IN  OUTPUT  EDIT 

I  PL  —1  -1-NO  PLOTS  REQUESTED 

1-CALL  FOR  OUTPUT  PLOTS 
ON  TAPE22 

ITPE-  0  -1-READ  FNK  FROM  TAPE24  AND 

COMPLETE  CALCULATION 
0-CREATE  FNK  ON  TAPE24  AND 
COMPLETE  CALCULATION 
1-CREATE  FNK  IN  TAPE24  AND  STOP 
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c)  Flux  Output:  File  01 


MINIMUM  CONTAINER  ARRAY  DIMENSION-  18324 

PLANE  GEOMETRY 

SOURCE: 

AT  X-0 
ISOTROPIC 

SCATTERING  ORDER  LL-  0 
LEGENDRE  COEFFICIENTS: 

1 .00000E+00 

ALA-  1 . OOOOE+OO 

BE  -  0 . 00000E+00  SI0-  1.00000E+00 
FLUX  CALCULATION 

S/X  1. 0000E+00  2. 0000E+00  3. OOOOE+OO  4. 0000E+00  5. 0000E+00 

♦★★★★★★A********************************************************************* 

1. OOOOE+OO  1 . 8394E-01  0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
2. OOOOE+OO  2 . 6256E-01  3.3834E-02  0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
3. OOOOE+OO  2 . 394 3E-01  9.3839E-02  8.2978E-03  0. OOOOE+OO  0. OOOOE+OO 
4. OOOOE+OO  2 . 1736E-01  1.1361E-01  3.3288E-02  2.2895E-03  0. OOOOE+OO 
5. OOOOE+OO  1.9957E-01  1.2105E-01  4.9595E-02  1.1823E-02  6.7379E-04 
6. OOOOE+OO  1.8523E-01  1.2322E-01  6.0605E-02  2.0710E-02  4.2106E-03 
7. OOOOE+OO  1.7348E-01  1.2293E-01  6.8028E-02  2.8447E-02  8.4158E-03 
8. OOOOE+OO  1. 6364 E- 01  1.2143E-01  7.3033E-02  3.4901E-02  1.2790E-02 
9. OOOOE+OO  1 . 5528E-01  1.1935E-01  7.6385E-02  4.0186E-02  1.7004E-02 
1 . 0000E+01  1.4806E-01  1.1700E-01  7.8590E-02  4.4479E-02  2.0905E-02 
1 . 1000E+01  1 . 4175E-01  1.1455E-01  7.9986E-02  4.7954E-02  2.4433E-02 
1 . 2000E+01  1 . 3617E-01  1.1209E-01  8.0803E-02  5.0760E-02  2.7583E-02 
1.3000E+01  1 . 3120E-01  1.0969E-01  8.1201E-02  5.3024E-02  3.0372E-02 
1.4000E+01  1 . 2674E-01  1.0737E-01  8.1293E-02  5.4845E-02  3.2828E-02 
1.5000E+01  1.2269E-01  1.0514E-01  8.1159E-02  5.6305E-02  3.4986E-02 
1.6000E+01  1 . 1901E-01  1.0300E-01  8.0860E-02  5.7470E-02  3.6876E-02 
1 . 7000E+01  1 . 1564E-01  1.0096E-01  8.043BE-02  5.8391E-02  3.8531E-02 
1 . 8000E+01  1 . 1254E-01  9.9019E-02  7.9927E-02  5.9111E-02  3.9978E-02 
1 . 9000E+01  1 . 0968E-01  9.7167E-02  7.9349E-02  5.9664E-02  4.1241E-02 
2 . 0000E+01  1.0702E-01  9.5402E-02  7.8725E-02  6.0078E-02  4.2344E-02 
***************************************************************************** 


NUMBER  OF  TERMS  FOR  CONVERGENCE 
00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO 
20 . 0180E-05  00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO 
18 . 0200E-05  20 . 0200E-05  00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO 
18 . 02 30E-05  18.0230E-05-22.0230E-05  00. OOOOE+OO  00. OOOOE+OO 
20 . 0260E-05  20.0260E-05  22 . 0260E-05-22 . 0260E-05  00. OOOOE+OO 
20 . 0280E-05  20.0280E-05  20 . 0280E-05-22 . 0280E-05-24 . 0300E-06 
20 . 0320E-05  22.0320E-05  22 . 0320E-05-24 . 0320E-05-24 . 0320E-05 
18 . 0310E-05  22.0330E-05  24.0330E-05  24 . 0330E-05-22 . 0330E-05 
20 . 0340E-05  24.0390E-05  24.0390E-05  22 . 0390E-05-24 . 0390E-05 
22 . 0360E-05  24.0370E-05  24.0370E-05  24 . 0370E-05-26 . 0380E-05 
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22 . 0380E-05 
22 . 0390E-05 
24 . 0410E-05 
24 . 0420E-05 
26 . 0450E-05 
26 . 0460E-05 
26 . 0470E-05 
28 . 0490E-05 
28 . 0510E-05 
28 . 0520E-05 


24 . 0390E-05 
26 . 0410E-05 
26 . 0420E-05 
26 . 04  30E-05 
26 . 0450E-05 
28 . 0470E-05 
28 . 0480E-05 
28 . 0490E-05 
28 . 0510E-05 
28 . 0520E-05 


22 . 0390E-05 
24 . 0410E-05 
26 . 0420E-05 
26 . 0430E-05 
26 . 0450E-05 
28 . 0470E-05 
28 . 0480E-05 
28 . 0490E-05 
30.0520E-05 
30 . 0530E-05 


26 . 0400E-05 
26 . 0410E-05 
26 . 0420E-05 
28.0440E-05 
28 . 0460E-05 
28 . 0470E-05 
26 . 0480E-05 
26 . 0490E-05 
28 . 0520E-05 
28 . 0530E-05 


26 . 0400E-05 
24 . 0410E-05 
26 . 0420E-05 
26 . 0440E-05 
28 . 0460E-05 
28 . 0470E-05 
30 . 0490E-05 
30 . 0500E-05 
30 . 0520E-05 
30 . 0530E-05 
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2)  Problem  2 


In  this  problem  a  beam  source  emitting  electrons  in  direction  p0  =  1  is  located  at  x  =  0.  A 
corrected  screened-Rutherford  (SR)  scattering  kernel  is  assumed  with  L  =  35  and  r?  =  0.05.  The 
points  as  in  problem  1  are  edited.  Note  the  appearance  of  the  unlikely  number  7.7777E+17  when 
This  value  indicates  the  delta  function  "infinity"  of  the  uncollided  flux. 


transport 

same 

s  =  x. 
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a)  Input:  File  SLEET.DAT 


SAMPLE  PROBLEM  2: BEAM  SOURCE  AT  X-O.SR-LL-35 
2  0  1.0  1.0 
35  -1  1 
0.05  1.0 
0.0  1.0 
20  5  1  1 
1 

90  90  3  3 
1.0E-04 
111-1 
0 

1.0  1.0 
1.0  1.0 
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b)  Flux  Output:  File  01 


MINIMUM  CONTAINER  ARRAY  DIMENSION-  26844 

PLANE  GEOMETRY 

SOURCE : 

AT  X-0 
ANISOTROPIC 
MU0-  1 . OOOOE+OO 

SCATTERING  ORDER  LL-  35 
LEGENDRE  COEFFICIENTS: 

1 . 00000E+00  8 . 59679E-01  7.07990E-01  5.68955E-01  4.49953E-01 

3.51820E-01  2 . 72761E-01  2.10075E-01  1.60941E-01  1.22766E-01 

9 . 33080E-02  7.07013E-02  5.34309E-02  4.02871E-02  3.03158E-02 

2.27721E-02  1.70787E-02  1.27907E-02  9.56708E-03  7.14768E-03 

5.33453E-03  3.97750E-03  2.96308E-03  2.20559E-03  1.64052E-03 

1 . 21938E-03  9 . 05761E-04  6.72402E-04  4.98886E-04  3.69953E-04 

2 . 74206E-04  2.03145E-04  1.50435E-04  1.11355E-04  8.23959E-05 

6 . 09456E-05 

ETA-  5 . 0000E-02  ALA-  1.0000E+00 
BE  -  0.00000E+00  SlO-  1.00000E+00 
FLUX  CALCULATION 

S/X  1 . 0000E+00  2 . 0000E+00  3.0000E+00  4.0000E+00  5.0000E+00 
***************************************************************************** 
1.0000E+00  7 . 7777E+17  O.OOOOE+OO  0.0000E+00  O.OOOOE+OO  0.0000E+00 
2 . OOOOE+OO  9.9542E-02  7.7777E+17  O.OOOOE+OO  0.0000E+00  0.0000E+00 
3 . OOOOE+OO  5 . 9803E-02  2.0083E-01  7.7777E+17  O.OOOOE+OO  O.OOOOE+OO 
4 . 0000E+00  4 . 9196E-02  1.0513E-01  2.8002E-01  7.7777E+17  O.OOOOE+OO 
5. 0000E+00  4 . 4716E-02  7.7885E-02  1.4534E-01  3.1745E-01  7.7777E+17 
6. OOOOE+OO  4 . 2405E-02  6.5544E-02  1.0361E-01  1.7253E-01  3.1605E-01 
7. 0000E+00  4 . 1050E-02  5.8607E-02  8.3991E-02  1.2251E-01  1.8462E-01 
8. OOOOE+OO  4 . 0167E-02  5.4178E-02  7.2725E-02  9.7921E-02  1.3318E-01 
9. OOOOE+OO  3 . 9526E-02  5.1074E-02  6.5419E-02  8.3464E-02  1.0649E-01 
1.0000E+01  3 . 9011E-02  4.8751E-02  6.0266E-02  7.3958E-02  9.0329E-02 
1 . 1000E+01  3 . 8560E-02  4.6911E-02  5.6409E-02  6.7214E-02  7.9516E-02 
1 . 2000E+01  3 . 8136E-02  4.5391E-02  5.3379E-02  6.2157E-02  7.1766E-02 
1 . 3000E+01  3 . 7725E-02  4.4087E-02  5.0914E-02  5.8200E-02  6.5927E-02 
1.4000E+01  3 . 7311E-02  4.2942E-02  4.8848E-02  5.4999E-02  6.1354E-02 
1 . 5000E+01  3 . 6894E-02  4.I913E-02  4.7076E-02  5.2344E-02  5.7654E-02 
1 . 6000E+01  3 . 6470E-02  4.0971E-02  4.5520E-02  5.0082E-02  5.4592E-02 
1 . 7000E+01  3.6043E-02  4.0100E-02  4.4143E-02  4.8128E-02  5.2004E-02 
1 . 8000E+01  3.5609E-02  3.9287E-02  4.2904E-02  4.6416E-02  4.9779E-02 
1.9000E+01  3 . 5174E-02  3.8521E-02  4.1776E-02  4.4897E-02  4.7839E-02 
2 . OOOOE+Ol  3 . 4737E-02  3.7797E-02  4.0742E-02  4.3533E-02  4.6129E-02 
***************************************************************************** 


NUMBER  OF  TERMS  FOR  CONVERGENCE 
00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO 
-53 . 0131E-05  00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO  00. OOOOE+OO 
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50.0161E-05-44.0161E-05 
46. 0181 E- 05- 40, 0181 E- 05 
39. 0201 E- 05- 37. 0201 E- 05 


-35 . 0221E-05 
-30 . 0251E-05 
80 . 0270E-05 
34 . 0290E-05 
34 . 0300E-05 
32 . 0330E-05 
30 . 0340E-05 
27 . 0390E-05 
29 . 0380E-05 
29 . 0390E-05 
30 . 0400E-05 
31. 0430 E- 05 
29 . 04  30E-05 
27 . 0450E-05 
27 . 0470E-05 


32 . 0221E-05 
29 . 0251E-05 
26 . 0270E-05 
34 . 0290E-05 
34 . 0300E-05 
31 . 0330E-05 
29 . 0340E-05 
28 . 0390E-05 
26 . 0380E-05 
30 . 0390E-05 
30 . 0400E-05 
30 . 0430E-05 
28 . 0430E-05 
28 . 0450E-05 
26 . 0470E-05 


00.0000E+00 
38 . 0181E-05 
35 . 0201E-05 
31 . 0221E-05 
86 . 0250E-05 
26 . 0270E-05 
33 . 0290E-05 
34 . 0300E-05 
31 . 0330E-05 
30 . 0340E-05 
28 . 0390E-05 
25.0380E-05 
27 . 0390E-05 
28.0400E-05 
27 . 0430E-05 
29 . 0430E-05 
26 . 0450E-05 
26.0470E-05 


00 . 0000E+00 
00 . 0000E+00 
34 . 0201E-05 
31 . 0221E-05 
77 . 0250E-05 
26 . 0270E-05 
33 . 0290E-05 
34 . 0300E-05 
31 . 0330E-05 
30 . 0340E-05 
27 . 0390E-05 
24 . 0380E-05 
23 . 0390E-05 
28 . 0400E-05 
30 . 0430E-05 
28 . 0430E-05 
27 . 0450E-05 
27 . 0470E-05 


00 . 0000E+00 
00 . 0000E+00 
00 . 0000E+00 
73 . 0220E-05 
64 . 0250E-05 
26 . 0270E-05 
33 . 0290E-05 
32 . 0300E-05 
30 . 0330E-05 
29. 0340 E- 05 
26 . 0390E-05 
25 . 0380E-05 
25 . 0390E-05 
29 . 0400E-05 
28 . 0430E-05 
28 . 0430E-05 
26 . 0450E-05 
26 . 0470E-05 
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3)  Problem  3 

The  final  problem  is  identical  to  problem  2  but  the  source  is  located  in  the  entire  left  half-space  (x 

<  0). 
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a)  Input:  File  SLEET.DAT 


SAMPLE  PROBLEM  3 : BEAM  SOURCE  7N  A  HALF-SPACE, SR-LL-3 5 
2  1  1.0  1.0 
35  -1  1 
0.05  1.0 
0.0  1.0 
5  22  1  1 
1 

75  75  3  3 
1.0E-04 
11-1-1 
0 

1.0  2.0 
-11.0  1.0 
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MINIMUM  CONTAINER  ARRAY  DIMENSION-  19018 


PLANE  GEOMETRY 


SOURCE: 

CONSTANT  HALF-SPACE 

ANISOTROPIC 

MU0-  1 . 0000E+00 


SCATTERING  ORDER  LL-  35 
LEGENDRE  COEFFICIENTS: 

1 .  OOOOOE-t-OO  8 . 59679E-01 
3 . 51820E-01  2 . 72761E-01 
9 . 3 3080E-02  7.07013E-02 
2 . 27721E-02  1.70787E-02 
5 . 33453E-03  3.97750E-03 
1 . 21938E-03  9 . 05761E-04 
2.74206E-04  2.03145E-04 
6 . 09456E-05 


7 . 07990E-01 
2 . 10075E-01 
5 . 34309E-02 
1 . 27907E-02 
2.96308E-03 
6.72402E-04 
1.50435E-04 


5 . 68955E-01 
1 . 60941E-01 
4 . 02871E-02 
9 . 56708E-03 
2 . 20559E-03 
4 . 98886E-04 
1 . 11 355E-04 


4 . 49953E-01 
1 . 22766E-01 
3.03158E-02 
7.14768E-03 
1.64052E-03 
3 . 69953E-04 
8 . 23959E-05 


ETA-  5 . 0000E-02  ALA-  1.0000E+00 


BE  -  0 . 00000E+00  SI0-  1.00000E+00 


FLUX  CALCULATION 


X/S  1.0000E+00  3 • OOOOE+OO  5.0000E+00  7. OOOOE+OO  9.0000E+00 
***************************************************************************** 


-1 . 1000E+01 
-1 . 0000E+01 
-9. 0000E+00 
-8. OOOOE+OO 
-7. OOOOE+OO 
6. 0000E+00 
-5-O000E+00 
-4. 0000E+00 
-3. 0000E+00 
-2. 0000E+00 
-1. 0000E+00 
0. OOOOE+OO 
1. OOOOE+OO 
2. OOOOE+OO 
3. OOOOE+OO 
4. OOOOE+OO 
5. OOOOE+OO 
6. OOOOE+OO 
7. OOOOE+OO 
8. OOOOE+OO 
9. OOOOE+OO 
1 . OOOOE+Ol 


1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
9.9354E-01 
3 . 6788E-01 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
1. OOOOE+OO 
9 . 9894E-01 
9 . 9352E-01 
9.7833E-01 
9 . 3998E-01 
8 . 2800E-01 
4 . 9787E-02 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


1. OOOOE+OO  1. OOOOE+OO  1. OOOOE+OO 
1. OOOOE+OO  1. OOOOE+OO  1. OOOOE+OO 
1. OOOOE+OO  1. OOOOE+OO  1. OOOOE+OO 
1. OOOOE+OO  1. OOOOE+OO  9.9995E-01 
1. OOOOE+OO  1. OOOOE+OO  9.9947E-01 
1. OOOOE+OO  9 . 9987E-01  9.9794E-01 
1. OOOOE+OO  9.9886E-01  9.9467E-01 
9 . 9962E-01  9.9600E-01  9.8884E-01 
9.9746E-01  9 . 9023E-01  9.7961E-01 
9 . 9184E-01  9.8029E-01  9.6591E-01 
9.8062E-01  9 . 6451E-01  9.4663E-01 
9 . 6043E-01  9.4087E-01  9.2044E-01 
9.2579E-01  9.0641E-01  8.8577E-01 
8.6627E-01  8.5712E-01  8.4066E-01 
7.5907E-01  7 . 8660E-01  7.8269E-01 
5 . 4290E-01  6 . 8482E-01  7.0860E-01 
6 . 7379E-03  5.3385E-01  6.1410E-01 
0. OOOOE+OO  3 . 0161E-01  4.9344E-01 
0. OOOOE+OO  9 . 1188E-04  3.3960E-01 
0. OOOOE+OO  0. OOOOE+OO  1.5027E-01 
0. OOOOE+OO  0. OOOOE+OO  1.2341E-04 
0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 


NUMBER  OF  TERMS  FOR  CONVERGENCE 
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OO.OOOOE+OO 
00 . 0000E+00 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
37 . 0100E-05 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
00 . OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
00 . OOOOE+OO 
00 . OOOOE+OO 


OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
00 . OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
23 . 0150E-05 
23 . 0150E-05 
23 . 0150E-05 
23 . 0150E-05 
24 . 0150E-05 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
00 . OOOOE+OO 
OO.OOOOE+OO 
00 . OOOOE+OO 
OO.OOOOE+OO 
00 . OOOOE+OO 


OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
22 . 0200E-05 
23 . 0200E-05 
24 . 0200E-05 
23 . 0200E-05 
23 . 0200E-05 
23 . 0200E-05 
24 . 0200E-05 
23 . 0200E-05 
26 . 0200E-05 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
00 . OOOOE+OO 
OO.OOOOE+OO 


OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
22 . 0260E-05 
22 . 0260E-05 
22 . 0260E-05 
22 . 0260E-05 
22 . 0260E-05 
21 . 0260E-05 
21 . 0260E-05 
21 . 0260E-05 
22 . 0260E-05 
22 . 0260E-05 
22 . 0260E-05 
23 . 0260E-05 
25 . 0260E-05 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 


OO.OOOOE+OO 
OO.OOOOE+OO 
OO.OOOOE+OO 
20 . 0290E-05 
20 . 0290E-05 
20 . 0290E-05 
19 . 0290E-05 
19 . 0290E-05 
18 . 0290E-05 
19 . 0290E-05 
20 . 0290E-05 
19 . 0290E-05 
20 . 0290E-05 
19 . 0290E-05 
20. 0290 E- 05 
20 . 0290E-05 
19 . 0290E-05 
21 . 0290E-05 
21 . 0290E-05 
21 . 0290E-05 
OO.OOOOE+OO 
OO.OOOOE+OO 
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C)  Programming  Notes 

The  SLEET  code  is  structured  on  a  modular  basis.  The  main  routine  calls  the  input  routines  and 
for  the  calculation  or  input  from  tape  (TAPE  24)  of  f^k  as  well  as  for  an  edit  of  the  input  options  (on 
TAPE  21  or  23)  if  desired.  In  addition,  the  path  is  advanced  and  a  call  for  the  flux  determination  (in 
SCALR)  is  made  along  with  the  output  edit.  The  use  of  dynamic  storage  allows  maximum  flexibility  in  the 
problem  size.  A  container  array  size  is  located  in  a  parameter  statement  in  the  main  program.  If  IASK.  > 
IMX  where  IMX  is  the  maximum  storage  requirement  for  the  given  problem,  a  diagnostic  message  is 
printed  on  the  screen  and  the  calculation  is  terminated.  The  minimum  container  array  required  is  also 
printed  in  the  output  (TAPE  21)  for  use  in  the  sizing  of  problems. 
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APPENDIX 

Derivation  of  Legendre  Scattering  Coefficients  for  the  Screened  Rutherford  Scattering  Kernel. 

The  derivation  begins  with  the  screened-Rutherford  scattering  law  modified  to  incorporate  a 
localized  electron-atom  collision  (screening  rj  >  0): 


«<«  )  .  -L  nLhnlD- 

8(/V  2*  (1+17  -  M0)2 


If  g  is  expanded  in  terms  of  Legendre  polynomials,  we  have 
oo 

1 

B) 


8(M(,)  =  X  wt  P*(M° 
1=0 


with 


f* 

=  2tr  j  dji0  P t(n0)  g(/i0) 
and  upon  substitution 


By  noting  that 
Ql(z)  = 


5  J.,d"” 


Z  =*  l+>? 


(A.l) 


(A.2) 


(A. 3) 


and 

dQi(z) 

dz 

Eq.  (A.3)  becomes 


1 1> 


(z-M0)2  ’ 


w t 


3(U3)^ 


dz 


Since 


d Qt  ^  Q\(z) 

dz  “  ’  S\-z> 

where  is  an  associated  Legendre  function  of  the  second  kind  of  order  t  satisfying 
1  Qii  *  (2<+lXi+*?)  -  (£+1)  Q}-i  , 


(A.4) 


(A. 5) 
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we  can  write  for  Eq.  (A.4) 


«*-  J  <!♦§) 


QtU) 

vTT= 


(A.6) 


From  Eq.  (A. 5)  therefore, 

iuj£+1  *  (2£+i)([+f?)  U£  -  (£+1)  W£_  1  .  (A. 7) 

This  recurrence  relation  is  initiated  by  w0,  w1  given  by  Eqs.  (5a,b)  which  are  obtained  directly  from  Eq. 

(A. 2). 
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PROGRAM  SLEET 

£**** 

C  MAIN  PROGRAM 
C  CHOOSES  INPUT  OPTIONS 

C  ADVANCES  PATH  COORDINATE 

C  CHOOSES  OUTPUT  EDIT  OPTIONS 

Qk  it  it  k 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

PARAMETER  (IASK-50000) 

COMMON/C1/IR , IS , JTI , LL , ITRN , LLP1 , LLP2 , ITC , IPM , LS , LX , JT, JX , 

*  IDFF , LN , LNP , LK , LKP ,  IEN , IEK , ITR ,  IP , IXY ,  IPL , ITPE ,  NS , 

*  KFF,KII ,LM, 1ST, IN, INF , NKT , I FLGl , NKX , IMXM , LAP , LOP , MLT 
COMMON/C2/IDEN( 18 ) , Q0 , AMO , ETA, ALA, BE , Bl , SIO , ERN , ERK , DC , SO , DS , X0 , DX 
COMMON/C 3/P I , X , S , Q , QY , TX , SXl , ET , TS , PHI , FD , DPH , WLPl 

COMMON/C 4 /I OTl , IOT2 , IOT3 , IOT4 , IOT06 , MX ( 2 ) , MY ( 2 ) , EXMX 
DIMENSION  A( 80000 ) 

OPEN ( UNIT- 20 , FILE- ' SLEETM . DAT ' , STATUS- ' OLD ' ) 

OPEN ( UNI T-21, FILE- 'Ol' , STATUS- ' UNKNOWN ' ) 

OPEN( UNIT-22,  FILE- '02' , STATUS- ' UNKNOWN ' ) 

OPEN( UN IT- 2 3 , FILE- '03 ' , STATUS- ' UNKNOWN ' ) 

OPEN( UNIT-24,  FILE- '04' ,  STATUS- ' UNKNOWN ' ) 

C ***  DIMENSION  OF  CONTAINER  ARRAY 
MLT- 3 

C*  *  * 

C  SET  CONSTANTS  AND  ARRAY  POINTERS 
C*CDC  EXMX- 199 
EXMX-80 
IN-20 
IOTl-21 
IOT2-22 
IOT3-23 
IOT4-24 
IOT06-6 
REWIND  IOT1 
REWIND  IOT2 
REWIND  IOT3 
REWIND  IOT4 
PI-4 . 0D0*DATAN( 1 . 0D0 ) 

IR-4 

WRITE ( IOT06, * )  'READ  INPUT  AND  SET  ARRAY  POINTERS' 

CALL  INPUT 
LA-LS*MX( 1 )+LX*MX( 2 ) 

LO-LS*MY( 1 )+LX*MY( 2) 

LNR-IR+LN 

LKR-IR+LK 

LLR-IR+LL 

LMR-LKR+MLT 

IXX-1 

ISS-IXX+LX 
IFTC-ISS+LS 
IAB-I FTC+LKR+MLT 
IOR-IAB+LA 
IDT-IOR+LO 
IMXT-IDT+LLR 
DO  1  1-1 , IMXT 
A( I ) -0 . 0 
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1  CONTINUE 

CALL  INPUTl ( A( IXX ) , A( I SS ) , A( I  FTC ) ,A( IAB ) ,A( IOR) 

*  ,  LX ,  LS  ,  LA,  LO  ,  LLR ,  LMR ) 

DO  2  I-IDT,IASK 

A( I )-0.0D0 

2  CONTINUE 
IGN-IDT 

IPO«IGN+LKR*LNR 

IFNT-IPO+LKR 

I F ( ITPE.LT.O)  GO  TO  103 

IHNT-IFNT+LNR*LKR 

IWT-IHNT+LNR*LKR 

IMX-IWT+LNR 

IMXl-IMX 

IF( IMX.LE.IASK)  GO  TO  150 
WRITE ( IOT06 ,  199  )  IASK, IMX 
WRITE ( IOT1 ,  199  )  IASK , IMX 

199  FORMAT( IX, ' INSUFFICIENT  STORAGE  IASK-' , 16 , 2X, 'NEED-' , 16 , '  FOR  COMP 
*LETION ' ) 

STOP 

150  CONTINUE 

WRITE ( IOT06 , * )  'DETERMINE  FLUX  EXPANSION  COEFFICIENTS' 

CALL  ALP4 ( A( I FNT ) ,A( IGN) ,A( IHNT) ,A( IWT) ,A( I FTC) ,A( IPO) 

*  , LNR , LKR ) 

WRITE ( IOT06 , * )  'EXPANSION  COEFFICIENTS  STORED  ON  FILE  02 ' 
IF(IP.LT.O)  GO  TO  112 
IF(IP.EQ.O)  IOT3-21 
I F  (  I P  .  EQ  .  1  ) 

♦WRITE ( IOT06 , * )  'INPUT  WRITTEN  ON  FILE  03' 

IF( IP.EQ.0) 

♦WRITE ( IOT06 , * )  'INPUT  WRITTEN  ON  FILE  01' 

CALL  INPRNT(A( IXX) , A( ISS) ,A( IFTC) ,LX,LS,LLR) 

IP— 1 

112  CONTINUE 

I F ( ITPE.NE. 0 )  STOP 
103  CONTINUE 

IF(IP.LT.O)  GO  TO  111 
IF(IP.EQ.O)  IOT3-21 
I F ( IP.EQ.l) 

♦WRITE ( IOT06 , ♦ )  'INPUT  WRITTEN  ON  FILE  03' 

IF(IP.EQ.O) 

*WRITE( IOT06 , * )  'INPUT  WRITTEN  ON  FILE  01' 

CALL  INPRNT(A( IXX) ,A( ISS ) , A( I  FTC ) , LX, LS , LLR ) 

111  CONTINUE 

DO  3  I-IGN, IASK 
A( I)-0.0D0 

3  CONTINUE 

I PHI -I FNT 

IERR-IPHI+LA*LO 

IQ-IERR+LA*LO 

IU-IQ+LKR 

IMX-IU+2*LKR 

I MX 2- I MX 

IMXM-IMX1 

I F ( I MX 2 . GT . IMX1 )  IMXM-IMX2 
IF( IMX. LE. IASK)  GO  TO  200 
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WRITE( IOTl , 199 )  IASK , I MX 
STOP 

200  CONTINUE 
LAP-LA 
LOP-LO 

DO  105  KST-1 , LS 
IST-KST 

S-A( ISS+KST-1 ) 

WRITE( IOT06 , * )  'BEGIN  S-',S 
SXl-JTI* ( S-l . 0D0 ) +1 . 0D0 
Q- S/ALA 

QY*( S/S 10 ) **BE 
QB-Q*QY/Bl 

IF(QB.LT.EXMX)  GO  TO  133 
LAP- ( KST-1 ) *MX { 1 ) +LX*MX ( 2 ) 

LOP-( KST-1 ) *MY ( 1 )+LX*MY( 2 ) 

GO  TO  104 
133  CONTINUE 

TS-DEXP ( -QB ) 

109  CONTINUE 

IF(S.EQ.O.O)  GO  TO  105 
TX-TS/S 

CALL  SCALR(A( IXX) ,A( IGN) ,A( IPHI ) ,A( I ERR ) ,A( IPO) ,A( IQ) 

*  ,  A( IFTC ) ,A( IU) , LNR , LKR , LA , LO , LX , LLR ) 

105  CONTINUE 

104  CONTINUE 

WRITE ( IOT06 , * )  'WRITE  OUTPUT  ON  FILE  01' 

CALL  OUT ( A( IAB ) ,A(IOR) ,A(IFTC) ,A( IPHI ) , A( I ERR ) , LA, LO , LLR ) 

STOP 

END 

SUBROUTINE  INPUT 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C  READ  AND  PREPARES  INPUT 

C*  *  *  * 

C  INPUT  ROUTINE 
FREE  FORMAT 

SUGGESTED  VALUES  ARE  GIVE  IN  PARENTHESIS 
**** 

COMMON/C 1/1 R , I S , JTI , LL , ITRN , LLPl , LLP2 , I TC , I PM , LS , LX , JT , JX , 

*  IDFF , LN , LNP , LK , LKP , I EN , IEK , ITR , I P , IXY , I PL , ITPE , NS , 

*  KFF,KII,LM,IST,IN, INF,NKT, IFLGl , NKX, IMXM, LAP , LOP , MLT 
COMMON/C2/I DEN ( 1 8 ) , Q0 , AMO , ETA , ALA , BE , B 1 , S I 0 , ERN , ERK , DC , SO , DS , X0 , DX 
COMMON/C  3/PI , X , S , Q , QY , TX , SX1 , ET , TS  f  PHI , FD , DPH , WLPl 

COMMON/C 4/1 OT1 , IOT2 , IOT3 , IOT4 , IOT06 , MX( 2 ) , MY( 2 ) , EXMX 

2#A#******A*#####AA*ir**#AA#**#*AA4rAA4r*AA#**#**AA***##*#**#A4rAAA** 

*  INPUT  DESCRIPTION 

-» 

:  CARD  1  (18A4) 

:  IDEN  -  PROBLEM  IDENTIFICATION 

:  CARD  2 

'  IS  -  1  ISOTROPIC  SOURCE 

-  2  ANISOTROPIC  SOURCE 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

w 

C 

c 

c 

c 

c 


JTI  -  0  DELTA  FUNCTION  SOURCE  AT 
X-0 

-  1  CONSTANT  SOURCE  IN  LEFT 

HALF-SPACE ( X< 0 ) 

Q0  -  SOURCE  STRENGTH 

AMO  -  COSINE  OF  ANGLE  OF  INCIDENCE  OF  SOURCE 
ELECTRONS 

CARD  3 

LL  -  ORDER  OF  SCATTERING  KERNEL  APPROXIMATION 

ITC  —1  SCREENED-RUTHERFORD  SCATTERING  KERNEL 

-  1  READ  IN  SCATTERING  COEFFICIENTS 

ITRN  —1  USE  TRANSPORT  CORRECTED  CROSS  SECTIONS 

-  1  DO  NOT  USE  TRANSPORT  CORRECTED  CROSS 

SECTIONS 

CARD  4 A  IF  ITC-1 
ALA  -  TOTAL  MEAN  FREE  PATH 


C  CARD  4B  IF  ITC— 1 

C  ETA  -  SCREENING  PARAMETER  FOR  SCREENED 
C  RUTHERFORD  SCATTERING  CROSS  SECTION 

C 

C  ALA  -  TOTAL  MEAN  FREE  PATH 


C 

C  CARD  5 

C  BE  -  POWER  OF  PATH  DEPENDENCE  OF  NUMBER 

C  OF  SECONDARIES  ( ALA( S )■( S/SIO ) **BE) 

C  SIO  -  CROSS  SECTION  PARAMETER 
C 

C  CARD  6 

C  LS  -  NUMBER  OF  PATH  GRID  POINTS 

C 

C  LX  -  NUMBER  OF  SPACE  GRID  POINTS 

C 

C  JT  —1 

-  1 

C 

C  JX  —1 

C  -  1 

C 

C  CARD  7 

IDFF— 1 

C 

C  -  1 

c 

C  CARD  8 

C  LN  -  MAXIMUM  NUMBER  OF  TERMS  IN 

MULTIPLE  COLLISION  EXPANSION 
C  OF  MOMENTS 

C 


READ 

IN 

PATH 

GRID 

POINTS 

READ 

IN 

PATH 

MESH 

SPACING 

READ 

IN 

SPACE 

GRID 

POINTS 

READ 

IN 

SPACE 

MESH 

SPACING 

COMPARISON  WITH  DIFFUSION  THEORY 
(BE-0  ONLY) 

SCALAR  FLUX  CALCULATION 
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C  LK  •  MAXIMUM  NUMBER  OF  TERMS  IN 

C  LEGENDRE  SERIES  EXPANSION 

C 

Z  IEK  -  NUMBER  OF  CONSECUTIVE  TERMS 
C  SATISFYING  RELATIVE  ERROR 

C  CRITERION  IN  LEGENDRE  EXP- 

C  SION ( 3 ) 

C  NOTE:IF  IEK.LT.0  SUPPRESS  ANALYTIC  DETERMINATION  OF 
C  OF  FIRST  COLLIDED  FLUX  WHEN  APPROPIATE 

r 

C  LM  -  NUMBER  OF  CYCLES  FOR  ACCELERATED 

C  CONVERGENCE ( 3 ) 

C 

C  CARD  9 

C  ERK  -  RELATIVE  TRUNCATION  ERROR  FOR 

C  LEGENDRE  EXPANSION  1 . 0E-04 ) 

C  CARD  10 

C  ITR  —1  PRINT  OUTPUT  ON  PRINTER  (132  COLUMNS) 

C  -  1  PRINT  OUTPUT  ON  TERMINAL  (72  COLUMNS) 

C 

C  IP  —1  DO  NOT  PRINT  INPUT  EDIT 

-  0  PRINT  INPUT  EDIT  WITH  OUTPUT  EDIT( TAPE21 ) 

C  -  1  PRINT  INPUT  EDIT  ON  TAPE23 

C 

C  IXY  —1  S  ON  ABCISSA, X  ON  ORDINATE  IN  EDIT 
C  -1  X  ON  ABCISSA, S  ON  ORDINATE  IN  EDIT 

C 

"  I  PL  —1  NO  PLOTS  REQUESTED 
C  -  1  WRITE  PLOT  FILE ( TAPE22 ) 

C 

C  CARD  1 1 

C  ITPE  —1  READ  FNK  FROM  TAPE24  AND  COMPLETE 
C  CALCULATION 

C.  rn  o  CREATE  FNK  ON  TAPE24  AND  COMPLETE 

CALCULATION 

C  -  1  CREATE  FNK  ON  TAPE24  AND  STOP 

C 

C  CARD  12  IF  ITC-1 

C  FTC  -  SCATTERING  COEFFICIENTS  <IC-0,LL) 

C 

CARD  13A  IF  JT— 1 

Z  SS(IT)-  PATH  GRID  POINTS ( IT-1 , LS ) 

Z  CARD  13B  IF  JT-1 

7  SO  -  INITIAL  PATH  LENGTH 

:  DS  -  PATH  INTERVAL 

:  CARD  14A  IF  JX— 1 

:  XX(IX)-  SPACE  GRID  POINTS  (IX»1,LX) 

:  CARD  14B  IF  JX-l 

:  X0  -  INITIAL  SPACE  GRID  POINT 

DX  -  SPACE  MESH  SPACING 

3**************************************************************** 
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C  READ  INPUT 

READ(IN,3)  IDEN 
3  FORMAT ( 18A4 ) 

READ( IN, * )  IS, JTI ,Q0,AM0 
READ ( IN , * )  LL , ITC , ITRN 
LLPl-LL+1 
LLP2-LL+2 

IF(ITC.LT.O)  READ ( IN , * )  ETA, ALA 

IF(ITC.GT.O)  READ (IN,*)  ALA 

READ( IN , *  )  BE, SIO 

I F ( BE . EQ .0.0)  SI0-1.0D0 

Bl-1 .ODO+BE 

KFF-2 

KII-2 

NS-1 

IF(IS.EQ.l)  GO  TO  106 

KFF-1 

KII-1 

IF(AMO.NE.O.O)  GO  TO  104 

KFF-2 

KII-2 

104  CONTINUE 
NS-2 

IF(JTI.GT.O)  NS-1 
IF(BE.NE.O.O)  NS-1 
106  CONTINUE 

READ( IN, * )  LS , LX, JT , JX 
READ( IN , * )  IDFF 
READ( IN, * )  LN , LK , IEK , LM 
IF(IEK.LT.O)  NS-1 
IEK-IABS( IEK) 

LNP-LN+1 

IF(LL.GT.LK)  LK-LL 

LKP-LK+1 

READ (IN,*)  ERK 

ERN-ERK/10.0D0 

IEN-IEK 

READ( IN, * )  ITR, IP, IXY, IPL 
READ( IN, * )  ITPE 
DO  70  J-1,2 
MY ( J ) -0 
MX( J )— 0 
70  CONTINUE 

IF( IXY.GT. 0 )  GO  TO  699 
MX( 1 ) -1 
MY( 2 )-l 
GO  TO  688 
699  CONTINUE 
MX ( 2 ) - 1 
MY ( 1 ) -1 
688  CONTINUE 
:  EXCLUSIONS 

IF(JTI.EQ.l)  LM-0 
IF(LL.EQ.O)  ITRN-1 
IF(BE.NE.O.ODO)  IDFF-1 
RETURN 
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END 

SUBROUTINE  INPUT1 ( XX , SS , FTC , AB , OR  f  LX , LS , LA , LO , LLR , LMR ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

READ  AND  PREPARE  FLOATING  POINT  INPUT 
*  *  * 

INPUT  ROUTINE 
FREE  FORMAT 

READ  AND  PREPARE  FLOATING  POINT  INPUT 
*  *  * 

COMMON/Cl/IR , IS , JTI , LL , ITRN , LLPl , LLP2 , I TC , I PM , LSS , LXX , JT , JX , 

*  IDFF , LN , LNP , LK , LKP , I EN , I EK , I TR , I P , IXY , I PL , I TPE , NS , 

*  KFF , KI I , LM , 1ST , IN , INF , NKT , I FLG1 , NKX , IMXM , LAP , LOP , MLT 
COMMON/C2/IDEN( 18 ) , Q0 , AMO , ETA, ALA, BE , Bl , SIO , ERN , ERK , DC , SO , DS , X0 , DX 
COMMON/C 3/P I , X , S , Q , QY , TX , SXl , ET , TS , PHI ,  FD , DPH , WLPl 

DIMENSION  XX ( LX ) , SS ( LS ) , FTC ( LLR ) , AB ( LA) ,OR(LO) 

WLP1-0.0D0 

IF(ITC.LT.O)  GO  TO  90 

READ( IN, * )  (FTC(IC) ,IC-IR,LLR) 

GO  TO  91 

90  CONTINUE 

CALL  SRUTH ( FTC , LMR ) 

91  CONTINUE 

I F ( ITRN. LT. 0.0)  CALL  TRNSC ( FTC , LLR ) 

F1-0.0D0 

IF(LL.GT.O)  Fl-FTC(IR+1) 

IF(Fl.NE.l.O)  DC— ALA/ ( 3 . 0D0* ( 1-Fl ) ) 

IF(JT.GT.O)  GO  TO  10 
READ( IN ,  *  )  ( SS ( IT ) , 1T«1 , LS ) 

GO  TO  20 

10  READ( IN, * )  SO , DS 

20  IF(JX.GT.O)  GO  TO  11 

READ( IN , * )  (XX(IX) ,IX-1,LX) 

GO  TO  21 

11  CONTINUE 

READ( IN, * )  X0 , DX 

21  IF(JT.LE.O)  GO  TO  60 
SS(1 )-S0 

IF(LS.EQ.l)  GO  TO  60 

DO  61  KST-2 , LS 

SS ( KST ) *SS ( KST-1 ) +DS 

61  CONTINUE 
60  CONTINUE 

IF(JX.LE.O)  GO  TO  62 
XX(1)-X0 

IF(LX.EQ.l)  GO  TO  62 
DO  63  IX-2,LX 
XX ( IX)*XX( IX-1 )+DX 
63  CONTINUE 

62  CONTINUE 

DO  50  IX-1, LX 
IF(IXY.LT.O)  OR ( IX ) -XX ( IX ) 

IF(IXY.GT.O)  AB( IX)-XX( IX) 

50  CONTINUE 

DO  51  IT-1, LS 
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IF(IXY.LT.O)  AB ( IT ) -SS ( IT ) 

IF(IXY.GT.O)  OR( IT)-SS( IT) 

51  CONTINUE 
RETURN 
END 

SUBROUTINE  INPRNT( XX , SS , FTC , LX , LS , LLR ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C*  *  *  * 

C  INPUT  PRINT  ROUTINE 
C  WRITES  INPUT  EITHER  ON  TAPE21 
C  OR  TAPE23 

Q*  *  *  * 

COMMON/Cl/IR , I S , JTI , LL , I TRN , LLPl , LLP2 , ITC , IFM , LSS , LXX , JT , JX , 

*  IDFF,LN,LNP,LK,LKP,IEN,IEK,ITR,IP,IXY,IPL,ITPE,NS, 

*  KFF , KI I , LM ,IST,IN,INF, NKT , I FLGl , NKX , IMXM , LAP , LOP , MLT 
COMMON/C2/IDEN ( 18) , Q0 , AMO , ETA, ALA, BE , Bl , SIO , ERN , ERK , DC , SO , DS ,X0 , DX 
COMMON/C 4/1 OT1 , IOT2 , IOT3 , IOT4 , IOT06 , MX ( 2 ) , MY ( 2 ) , EXMX 

DIMENSION  XX ( LX ) , SS ( LS ) , FTC ( LLR ) 

3  FORMAT ( 18A4 ) 

WRITE( IOT3 , 3 )  IDEN 
WRITE( IOT3 , 103  )  IS 

103  FORMAT ( /6X ,  '  I S  «', 12 , 14X, ' 1-ISOTROPIC  SOURCE'/ 

*  6X,20X, '2-BEAM  SOURCE') 

I F ( IS . EQ . 2 )  WRITE( IOT3 ,  121 )  AMO 

121  FORMAT(/5X, 'AMO  -', 1PE12 . 4 , 6X ,  ' COSINE  OF  ANGLE  OF  INCIDENCE'/ 

+  28X, 'SOURCE  PARTICLES') 

WR I TE ( I OT  3,200)  JTI 

200  FORMAT (/6X, 'JTI- ' ,12,1 4X, ' 0 -DELTA  FUNCTION  SOURCE  AT  X-0'/ 

+  26X, '1-CONSTANT  SOURCE  IN  LEFT  HALF-SPACE (X<0 ) ' ) 

WRITE ( IOT3 ,  201 )  Q0 

201  FORMAT (/6X, 'Q0  -', 1PE12 . 4 , 6X ,' SOURCE  STRENGTH') 

IF(ITC.GT.O)  GO  TO  199 

WRITE ( IOT3 , 204 )  ETA, ALA 

204  FORMAT (/5X, 'ETA  1PE12 . 4 , 6X ,' SCREENING  PARAMETER'/ 

*  5X, '  ALA-' ,1PE12.4,6X, 'MEAN  FREE  PATH') 

GO  TO  198 

199  CONTINUE 

WRITE ( IOT3 , 197 )  ALA 

197  FORMAT(/5X, 'ALA  - ' , 1PE12 . 4 , 6X, ' MEAN  FREE  PATH') 

198  CONTINUE 

WRITE( IOT3 , 202 )  LL 

202  FORMAT ( /6X , ' LL  - ' , 1 2 , 16X, ' ORDER  OF  SCATTERING  APPROXIMATION') 

WRITE ( IOT3 , 208 )  ITC,ITRN 

208  FORMAT (/5X, 'ITC  -',I2,12X,'  -1-USE  TEST  KERNEL'/ 

+  25X, '  1-READ  COEFFICIENTS'/ 

+  5X, ' ITRN- ' , 12 , 12X, '  -1-USE  TRANSPORT  CORRECTED  CROSS  SECTIONS'/ 

+  25X, '  1-DO  NOT  USE  TRANSPORT  CORRECTION') 

WRITE( IOT3 , 203 )  LL, ( FTC( IC) , IC»IR,LLR) 

203  FORMAT (/6X, 'SCATTERING  COEFFICIENTS ( L-0-' , 12 ,') '/ 

+  ( IX, 5( 1PE12 . 4 ) ) ) 

WRITE( IOT3 , 207 )  BE, SIO 

207  FORMAT ( /6X , ' BE  -', 1PE12 . 4 , 6X, ' EXPONENT  IN  VARIATION  OF  C'/ 

+  28X, 20H( C( S )-C* ( S/S 10 ) **BE )/ 

+  5X,'SI0  -' ,1PE12. 4, 6X, 'PARAMETER  IN  VARIABLE  C' ) 

WRITE ( IOT3 , 106 )  LS,LX 

106  FORMAT ( /6X , ' LS  - ' , 1 3 , 1 5X, ' NUMBER  OF  PATH  POINTS'/ 


A45 


May  30  11:35  1990  sleetm.for  Page  9 


*  6X, ' LX  -' ,I3,15X, 'NUMBER  OF  SPACE  POINTS') 

WRITE ( IOT3 , 105 )  JT,JX 

105  FORMAT ( /6X, 'JT  -',I2,12X,'  -1-READ  IN  PATH  POINTS'/ 

*  24X, '  1-READ  IN  MESH  SPACING'/ 

*  6X,'JX  -' , 12 , 12X, '  -1-READ  IN  SPACE  POINTS'/ 

*  24X, '  1-READ  IN  MESH  SPACING' ) 

ZP( JT.GT.O) 

+  WRITE ( IOT3 , 117 )  S0,DS 

117  FORMAT {/6X, 'SO  -', 1PE12 . 4 , 6X, ' INITIAL  PATH'/ 

*  6X,'DS  -' ,1PE12. 4,6X, 'PATH  INTERVAL') 

I F ( JX.GT.O) 

+  WRITE ( IOT3 , 118 )  X0,DX 

118  FORMAT ( 6X , ' X0  - ' , 1PE12 . 4 , 6X , ' FIRST  POSITION'/ 

*  6X,'DX  -' ,1PE12.4,6X, 'SPACE  INTERVAL') 

WRITE( IOT3 , 107 )  ( SS ( IT ) , IT-1 , LS ) 

107  FORMAT(/6X, 'PATH  POINTS '/( IX , 5 (1PE12 . 4 ) ) ) 

WRITE ( IOT3,108)  (XX(IX) ,IX-1,LX) 

108  FORMAT(/6X, 'SPACE  POINTS '/( IX, 5{  1PE12 . 4 ) ) ) 

WRITE( IOT3 , 300 )  IDFF 

300  FORMAT (/5X, ' IDFF-' , 12 , 12X, '  -1-DIFFUSION  COMPARISON'/ 

+  25X, '  1-SCALAR  FLUX  CALCULATION') 

WRITE( IOT3, 303)  LN 

303  FORMAT ( /6X, ' LN  - ' , 1 3 , 15X , ' MAXIMUM  NUMBER  OF  TERMS  IN'/ 

+  28X, 'MULTIPLE  COLLISION  EXPANSION') 

WRITE ( IOT3 , 304 )  ERN , IEN 

304  FORMAT ( 5X, 'ERN  -', 1PE12 . 4 , 6X, ' RELATIVE  TRUNCATION  ERROR'/ 

+  28X , ' FOR  MULTIPLE  COLLISION  EXPANSION'/ 

+  5X, ' IEN  -' ,12, 16X, 'NUMBER  OF  TERMS  SATISFYING  RELATIVE'/ 

+28X, 'ERROR  CRITERION') 

WRITE ( IOT3 , 110 )  LK 

110  FORMAT ( /6X , ' LK  - '  , 1 3 , 15X, ' MAXIMUM  NUMBER  OF  TERMS  IN'/ 

*  28X, 'LEGENDRE  EXPANSION') 

WRITE ( IOT3 , 111 )  ERK , IEK , LM 

111  FORMAT ( 5X, 'ERK  -', 1PE12 . 4 ,6X, 'RELATIVE  TRUNCATION  ERROR'/ 

*  5X,'IEK  -' ,12, 16X, 'NUMBER  OF  TERMS  SATISFYING  RELATIVE'/ 

*  28X, 'ERROR  CRITERION'/ 

*  5X, '  LM  -' ,12,1 6X, 'NUMBER  OF  PATTES  FOR  CONVERGENCE ( JTI-0 ) 

*'  ) 

WRITE ( IOT3 , 116 )  ITR,IP 

116  FORMAT (/5X, 'I TR  «',I2,12X,'  -1-OUTPUT  TO  PRINTER(132  COLUMNS)'/ 

*  26X, '1-OUTPUT  TO  TERMINAL ( 7 2  COLUMNS)'/ 

*/6X, ' IP  -' , 12 , 12X, '  -1-DO  NOT  PRINT  INPUT  EDIT'/ 

*  25X, '  0-PRINT  INPUT  EDIT  ON  TAPE21 '/ 

*  26X, '1-PRINT  INPUT  EDIT  ON  TAPE23 ' ) 

WRITE ( IOT3 , 305 )  IXY 

305  FORMAT (/5X, 'IXY  -',I2,12X,'  -1-S  ON  ABCISSA  AND  X  ON'/ 

+28X,  'ORDINATE  IN  OUTPUT  EDIT'/ 

+  25X, '  1-X  ON  ABCISSA  AND  S  ON'/ 

+  28X, 'ORDINATE  IN  OUTPUT  EDIT') 

WRITE( IOT3 , 307 )  IPL 

307  FORMAT (/5X, 'IPL  -',I2,12X,'  -1-NO  PLOTS  REQUESTED'/ 

*  25X, '  1-CALL  FOR  OUTPUT  PLOTS'/ 

*  28X, 'ON  TAPE22 ' ) 

WRITE( IOT3 , 306 )  ITPE 

306  FORMAT ( /5X, ' ITPE- ' , 12 , 12X , '  -1-READ  FNK  FROM  TAPE24  AND'/ 

+  28X, 'COMPLETE  CALCULATION'/ 
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+  25X,  '  0-CREATE  FNK  ON  TAPE24  AND  '/ 

+  28X, 'COMPLETE  CALCULATION'/ 

+  25X,  '  1 "CREATE  FNK  IN  TAPE24  AND  STOP') 

RETURN 

END 

SUBROUTINE  ALP4 ( FN,GN,HN,W, FTC , P , LNR , LKR ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/C 1/IR, IS, JTI , LL , ITRN , LLP1 , LLP2 , ITC, IPM, LS , LX, JT, JX, 

*  IDFF , LN , LNP , LK , LKP , IEN , IEK , ITR , IP , IXY , IPL , ITPE , NS , 

*  KFF , KI I , LM , 1ST , IN , INF , NKT, I FLG1 , NKX , IMXM , LAP , LOP , MLT 
COMMON/C2/IDEN ( 18 ) , Q0 , AMO , ETA, ALA, BE , B1 , SI  0 , ERN, ERK , DC, SO , DS , X0 , DX 
COMMON/C 4/1 OT1 , IOT2 , IOT3 , IOT4 , IOT06 , MX ( 2 ) , MY ( 2 ) , EXMX 

DIMENSION  GN ( LNR , LKR ) , FN( LNR , LKR ) , HN ( LNR , LKR ) , W( LNR ) , FTC ( LKR ) 

*  , P ( LKR ) 

FN(IR,IR)-1.0D0 

W(IR)-1.0D0 

I F( IS . EQ . 2 )  GO  TO  40 
FNl-1 .0D0 
Xl-FTC( IR )/Bl 
DO  20  N-1,LN 
NR-IR+N 
FNl»Xl*FNl 
FN ( NR , IR ) -FNl 
W<NR)-FN(NR, IR) 

20  CONTINUE 

GN( IR, IR+1 )»1 . 0D0/3 . 0D0 
GO  TO  62 

40  CONTINUE 

CALL  LEG (P, AMO, LK, LKR) 

DO  59  LXX-1, LKP 

L-LXX-1 

LR-IR+L 

FN ( IR.LR)-P(LR) 

59  CONTINUE 

W  ( I R )  -  FN  (  I R ,  I R ) 

DO  60  N-1,LN 
NR-IR+N 

DO  61  LXX-1 , LLPl 

L-LXX-1 

LR-IR+L 

Xi-FTC( LR)/Bl 

FN( NR , LR ) • ( Xl * *N) *P( LR) 

61  CONTINUE 

W ( NR ) - FN ( NR , I R ) 

60  CONTINUE 

41  CONTINUE 

DO  71  LXX-1, LKP 

L-LXX-1 

LR-IR+L 

GN ( I R , LR ) -AMO  *  P { LR ) 

71  CONTINUE 

62  CONTINUE 

DO  30  N-l , LN 
NR-IR+N 

DO  31  LXX-1, LKP 
L-LXX-1 
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LR-IR+L 

Tl-  ( L+l .  0D0  )  *FN(  NR ..  LR+1 ) 

T2-L*FN(NR,LR-1 ) 

T3-N*FTC ( LR  >  *GN( NR-1 , LR ) 

GN(NR,LR)-( (Tl+T2)/(2*L+1.0D0)+T3)/(N*Bl+1.0D0) 
31  CONTINUE 

W(NR)-GN(NR, IR) 

30  CONTINUE 

WRITE( IOT4 , * )  ( W ( NX ) , NX* I R , LNR ) 

DO  100  K-2.LK 
KR-IR+K 
AK2-2+K-1 
LU-LKP-K 

IF(IS.EQ.l)  GO  TO  52 
DO  51  LXX-1,LKP 
L-LXX-1 
LR-IR+L 

HN( IR,LR)-P(KR)*P(LR) 

51  CONTINUE 
GO  TO  53 

52  CONTINUE 

DO  54  LXX-1,LKP 
L-LXX-1 

HN( IR, IR+L ) -0 . 0D0 
54  CONTINUE 

HN( IR, KR)-1 . 0D0/( AK2+2 . 0D0 ) 

53  CONTINUE 

DO  101  N— 1 , LN 
NR-IR+N 

DO  102  LXX— 1 , LU 
L-LXX-1 
LR-IR+L 
L2-2*L+1 

Tl- { N*Bl+l . 0D0-K ) *FN{ NR, LR) 

T2-L*GN ( NR, LR-1 ) 

T3- ( L+l . 0D0 ) *GN( NR , LR+1 ) 

SS3- ( T2+T3 ) +AK2/L2 

T4-N*  FTC ( LR ) * ( HN ( NR-1 , LR ) -FN ( NR- 1 , LR ) ) 

SS5-T1+SS3+T4 

HN(NR,LR)-SS5/(K+N*B1 ) 

102  CONTINUE 

W ( NR ) -HN ( NR , I R ) 

101  CONTINUE 

WRITE( IOT4 , * )  (W(NX) , NX-IR,LNR) 

DO  201  NX— 1 ( LNP 
N— NX-1 
NR-IR+N 

DO  202  LXX-1 , LU 

L-LXX-1 

LR-IR+L 

FN(NR,LR)-GN(NR,LR) 

GN ( NR , LR ) -HN ( NR , LR ) 

202  CONTINUE 
201  CONTINUE 
100  CONTINUE 
RETURN 
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END 

SUBROUTINE  SCALR( XX, W, PHI , AMERR , PE , QU , FTC ,U , LNR , LKR, LA, LO , LX , LLR ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C*  *  *  * 

C  NUMERICAL  ALGORITHM 

C  CALCULATES  FLUX  DISTRIBUTION  PHI 

C  FOR  ALL  OPTIONS 

C  PREPARES  ARRAYS  FOR  OUTPUT  EDIT 

C  *  *  *  * 

COMMON/C 1/1 R , I S , JTI , LL , I TRN , LLPl , LLP2 , I TC , I PM , LS , LXX , JT , JX , 

*  IDFF , LN , LNP , LK , LKP , I EN , I EK , I TR , I P , IXY , I PL, ITPE,NS, 

*  KFF,KII,LM,IST,IN,INF,NKT,IFLG1, NKX , IMXM , LAP , LOP , MLT 
COMMON/C2/I DEN ( 18 ) , Q0 , AMO , ETA , ALA , BE , Bl , S 1 0 , ERN , ERK , DC , SO , DS , X0 , DX 
COMMON/C  3/P I , X , S , Q , QY , TX , SXl , ET , TS , PHI , FD , DPH , WLPl 

COMMON/C 4/1 OTl , IOT2 , IOT3 , IOT4 , IOT06 ,MX( 2 ) , MY ( 2 ) , EXMX 
COMMON/C5/  PHO , ETS , ETl , A2T 

DIMENSION  W ( LNR , LKR ) ,PHI(LA,LO) , AMERR ( LA , LO ) ,XX(LX) , FE ( LKR ) 

*  , QU ( LKR ) , FTC ( LLR ) , U ( 2 , LKR ) 

DATA  KFLG/1/ 

IF(KFLG.LT.O)  GO  TO  181 
REWIND  IOT4 

DO  180  K-KFF , LK , KI I 
KX-IR+K 

I F ( KI I . EQ . 2 )  READ ( IOT4 , * )  (  W( NX , 1 ) , NX-IR , LNR ) 

READ( IOT4 , * )  ( W ( NX , KX ) , NX— I R , LNR ) 

180  CONTINUE 
KFLG— 1 

181  CONTINUE 
ERNT-ERN 
NKX-0 

AMXN-1 . 0E-03 
DO  1001  IX-1,LX 
NDI--DLOG10 ( ERK ) 

IMXX-1 
SKL-1 .0D10 
X-XX( IX) 

WRITE ( IOT06 , * )  '  BEGIN  X-',X 

IA-MX( 1 ) *IST+MX( 2 ) *  IX 

IO-MY( 1 ) *IST+MY( 2 ) *IX 

ET-X/S 

ETS-1-ET*ET 

ETl-l-ET 

DET-DABS( ET) 

I F ( DABS ( DET-1 . 0D0 ) . LT . 1 . 0E-1 2 )  GO  TO  976 
IF(DET.GT.l.ODO)  GO  TO  977 
INF-1 

CALL  ANACOL ( FTC , QU , PE , LKR , LLR ) 

IF(INF.LT.O)  GO  TO  979 
AU-DEXP { QY*Q/Bl ) -1 . 0D0-FD 
W( 1 , IR ) -AU 
U ( 1 , I R ) -AU 

CALL  LEG( PE, ET, LK, LKR) 

155  CONTINUE 
AERS-ERN 
A2T-1 .0D10 
A2-A2T 
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ALT-A2T 

ICE-0 

SK-PH0+PH1+AU/2 . 0D0 
IFLGl-1 

IF(IMXX.LT.O)  NKX-0 
M-0 

DO  140  K-KFF , LK ,  KI I 
KX-IR+K 

IF(Q.EQ.O.ODO)  GO  TO  172 

IF(K.LE.NKX)  GO  TO  152 

A1-1.0D10 

SN-0 . 0D0 

IE-0 

TZ-1.0D0 

IF(NS.EQ.2)  TZ«Q*QY 
DO  150  N-NS , LN 
TZ«Q*QY*TZ/N 
TN-TZ*W( IR+N , KX ) 

SN-SN+TN 

C*** 

C  DIAGNOSTIC  PRINT  1 
C  WRITE (21,97)  N , TN , SN , Al 

97  FORMAT ( 1X,I4,5(1PE12.4) ) 

IF(SN.NE.O.O)  Al-DABS ( TN/SN ) 

IF(A1 .LT.ERN)  IE-IE+1 

AERN-N* 1 . OE-3 

IF( IE.EQ.IEN)  GO  TO  151 

150  CONTINUE 

151  CONTINUE 
W( 1,KX)-SN 

152  CONTINUE 

IF(AERN.GT.AMXN)  AMXN-AERN 
SR- ( W( 1 , KX ) -W ( 1 , KX-2 ) * JTI ) 

PKK«PE( KX-JTI ) 

TK-( 2*K* ( 1- JTI ) +1 . 0D0 ) *SR*PKK/2 . 0D0 

W(  2 , KX ) — TK 

U( 1 , KX ) -SR 

SK-SK+TK 

W( 3,KX)-SK 

IF(SK.NE.O.ODO)  A2-DABS ( TK/SK ) 

IF ( A2 . LT . ERK . AND . ALT . LT . ERK )  IEE-IEE+1 
ALT-A2 

I F ( A2 . GT . A2T . OR . K . LT . LK/4 )  GO  TO  142 
NKT-K 

I F ( A2 . GT . 1 . OE-1 5 )  A2T-A2 
142  CONTINUE 
•  ■»  *  * 

:  DIAGNOSTIC  PRINT  2 
:  WRITE( 21,97)  K , TK , SK , A2 

'*** 

NKX-K 

IF( IEE.EQ. IEK)  GO  TO  141 
140  CONTINUE 

IF(SK.EQ.O.O)  GO  TO  172 
IFLG1— 1 
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SK-W<  3, IR+NKT) 

141  CONTINUE 

TlL*DABS ( W( 2 , IR+KFF ) ) 

DO  161  K-KFF,NKX,KII 
Tl-DABS{W( 2, IR+K) ) 

IF(TI.GT.TIL)  TlL-Tl 

161  CONTINUE 
TlL-TlL/DABS ( SK ) 

IF(TlL.LE.l.ODO.OR.TlL.EQ.O.O)  GO  TO  162 
ATI-DLOGIO(TIL) 

El-DABS ( (SK-SKL)/SK) 

IF(El.LT.ERK)  GO  TO  162 

SKL-SK 

IMXX— 1 

LMl-ATl+l+NDI 

LM2— DLOGIO(ERN) 

IF ( LM2 . GT . LMl . OR . LMl . GT . 20 )  GO  TO  162 
ERN-1 . ODO/IO . 0D0**LM1 
NDI-NDI+1 
GO  TO  155 

162  CONTINUE 

I F ( I FLG1 . GT . 0 . OR . LM . EQ . 0 )  GO  TO  163 
CALL  ACC ( U , W , PE , SK , LNR , LKR , M ) 

163  CONTINUE 
ERN-ERNT 

IF(IFLGl.GT.O)  GO  TO  172 
XI— DLOG10(A2T)+1 
IX1-X1 

AERS-10.0**IXl 
172  CONTINUE 

IF(IFLGl.LT.O)  NKX-NKT 
IF(NKX.LT.IO)  NKX-10 

AMERR( IA , 10 ) ■ ( NKX+AMXN+M* 1 . 0E-04 ) *IMXX*AERS 
PHI ( IA, IO)-TX*SK*SX1*QO 
SKP-PHI ( IA, 10) 

IF( IDFF.GT.O)  GO  TO  1001 
CALL  DIFF 

PHI ( IA, 10) -DABS ( ( SKP-DPH ) /SKP ) 

135  CONTINUE 
GO  TO  1001 
976  CONTINUE 

IF(JTI.EQ.l)  GO  TO  87 
PHI( IA, IO)-TX*Q0/2.0D0 
IF(IS.EQ.l)  GO  TO  978 
PHI( IA, IO)-O.ODO 

I F ( DABS( ET*AM0-1 . 0D0 ) . LT. 1 . OE-12 )  GO  TO  979 
GO  TO  978 
87  CONTINUE 

PHI( IA,IO)-O.ODO 

I F  (  I S  .  EQ .  2  .  AND  .AM0.EQ.1.0D0)  PHI«<IA,IO)-Q0*TS 
IF(ET.LT.O.O)  PHI ( IA, IO)-QO 
978  CONTINUE 

IF( IDFF.GT.O)  GO  TO  1001 
IF(PHI(IA,IO). EQ .0.0)  GO  TO  1001 
CALL  DIFF 

PHI ( IA, 10) -DABS ( ( PHI ( IA, 10 ) -DPH )/PHI ( IA, 10) ) 
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GO  TO  1001 
977  CONTINUE 

IF(ET.GT.1.0.OR. JTI.EQ.O)  GO  TO  1001 
PHI { IA, IO)-Q0 
GO  TO  978 
979  CONTINUE 

PHI( IA, IO)«7 .7777E17 
IF( IDFF.GT.0)  GO  TO  1001 
PHI (IA,IO)“1.0D0 
1001  CONTINUE 
RETURN 
END 

SUBROUTINE  ANACOL ( FTC , QU , PE , LKR , LLR ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/C1/IR ,  IS , JTI , LL , ITRN , LLP1 , LLP2 , I TC , I PM , LS , LXX , JT , JX , 

*  IDFF , LN , LNP , LK , LKP , IEN ,  IEK , ITR ,  IP , IXY ,  I  PL , ITPE , NS , 

*  KFF , KI I , LM , I ST , IN , INF , NKT , I FLG1 ( NKX , IMXM , LAP , LOP , MLT 
COMMON/C2/IDEN( 18) , Q0 ,  AMO , ETA, ALA, BE , B1 , SI0 , ERN, ERK , DC , SO , DS , X0 , DX 
COMMON/C  3/P I ,X,S,Q,QY,TX, SX1 , ET , TS , PHI , FD , DPH , WLPl 

COMMON/C 4/1 OTl , IOT2 , IOT3 , IOT4 , IOT06 , MX( 2 ) , MY( 2 ) , EXMX 
COMMON/C 5/  PH0,ETS,ET1,A2T 
DIMENSION  FTC(LLR) ,QU(LLR) , PE(LKR) 

IF(JTI.EQ.l)  GO  TO  10 
IF(IS.EQ.l)  GO  TO  11 
PH0-0.0D0 

IF(ET.NE.AMO)  GO  TO  200 
INF— 1 
GO  TO  200 
11  CONTINUE 
PH0-0.5D0 
GO  TO  200 
10  CONTINUE 

IF(IS.EQ.l)  GO  TO  15 
PH0-0.0D0 

IF(AMO.GE.ET)  PH0-1.0D0 
GO  TO  200 
15  CONTINUE 

PH0-ET1/2.0D0 
200  CONTINUE 
PH1-0.OD0 
FD-0.0D0 

I F( NS . NE . 2 )  RETURN 

CALL  FCOLA( FTC, QU, PE , LKR, LLR ) 

FD-Q 

RETURN 

END 

SUBROUTINE  FCOLA( FTC , QU , PE , LKR , LLR ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

'*  *  *  * 

:  FIRST  COLLIDED  FLUX  FOR  IS-2 
DETERMINES  FIRST  COLLIDED  FLUX 
FOR  BEAM  SOURCE 
**** 

COMMON/Cl/IR,IS, JTI ,LL, ITRN , LLPl , LLP2 , ITC , IPM, LS , LX, JT, JX, 

*  IDFF , LN , LNP , LK , LKP , IEN, IEK , ITR, IP , IXY , IPL, ITPE , NS , 

*  KFF, KI I, LM, 1ST, IN, INF, NKT, I FLG1, NKX, IMXM, LAP, LOP, MLT 
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COMMON/C2/I DEN (18) , Q0 , AMO , ETA, ALA, BE, B1 , SI  0 , ERN , ERK , DC , SO, DS ,X0 , DX 
COMMON/C  3/P I , X , S , Q , QY , TX , SXl , ET , TS , PHl , FD , DPH , WLPl 
COMMON/C 4 /I OTl , IOT2 , IOT3 , IOT4 , IOT06 , MX( 2 ) , MY( 2 ) , EXMX 
DIMENSION  FTC(LLR) ,QU(LLR) ,PE(LKR) 

CALL  LEG( PE , AMO , LL , LKR ) 

DO  2  KP-l,LLPl 
K-KP-1 

QU(  IR+K)«PE( IR+K) 

2  CONTINUE 

CALL  LEG(PE,ET,LL, LKR ) 

IMF-1 

IF(AMO.LT.ET)  IMF-0 
IMS-1 

IF(AMO.GT.ET)  IMS-0 
Dl-DABS ( ET-AMO ) 

IF( Dl . LT. 1 . OE-12 )  GO  TO  101 

D2-1 .0D0-AM0 

T1-0.0D0 

T2-0.0D0 

I F ( D2 . GT . 0 . 0D0 )  Tl-DLOG( D2/D1 ) 

D3-1 .0D0+AM0 

I F ( D3 . GT . 0 . 0D0 )  T2-DLOG( D3/D1 ) 

P0-IMF*T2+IMS*T1 

Pl-AM0*P0- ( 1 . ODO  +  ET ) * IMF+ ( 1 .ODO-ET) *  IMS 

SL-FTC ( IR ) *P0/2 . 0D0 

IF(LL.EQ.O)  GO  TO  100 

SL-SL+3 . 0D0*FTC( IR+1 ) *AM0*Pl/2 . 0D0 

IF(LL.EQ.l)  GO  TO  100 

DO  1  L-2 , LL 

Tl-( 2*L-1 )*AM0*P1/L 

T2- ( L-l ) *P0/L 

T3-( PE( IR+L ) —PE ( IR+L-2 ) )/L 

P2-T1-T2-T3 

P0-P1 

P1-P2 

SL-SL+(2*L+1)*FTC( IR+L ) *QU ( IR+L ) *P2/2 . 0D0 
1  CONTINUE 

100  CONTINUE 
PHl-Q*SL 
RETURN 

101  CONTINUE 
INF— 1 
RETURN 
END 

SUBROUTINE  ACC ( U , W, PE , SK , LNR , LKR , MU ) 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z) 

COMMON/Cl/IR , I S , JTI , LL , ITRN , LLPl , LLP2 , ITC , I PM , LS , LXX , JT , JX , 

*  IDFF, LN, LNP, LK, LKP, IEN, IEK, ITR, IP, IXY, IPL, ITPE,NS, 

*  KFF , KI I , LM , 1ST , IN , INF , NKT , I FLGl , NKX , IMXM , LAP , LOP , MLT 
COMMON/C2/I DEN ( 1 8 ) , QO , AMO , ETA , ALA , BE , Bl , S 1 0 , ERN , ERK , DC , SO , DS , XO , DX 
COMMON/C 3/P I , X , S , Q , QY , TX , SXl , ET , TS , PHI , FD , DPH , WLPl 

COMMON/C 4/1 OTl , IOT2 , IOT3 , IOT4 , IOT06 , MX( 2 ) , MY( 2 ) , EXMX 

COMMON/C5/  PHO , ETS , ETl , A2T 

DIMENSION  W ( LNR , LKR ) , U( 2 , LKR ) , PE ( LKR ) 

A2ML-1 .0D10 
NKTL-NKT 
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IFLG1-1 

A2T-1 .0D10 

ALT-A2T 

LKMM-LK 

IEE-0 

SKL-SK 

DO  500  M«1,LM 
MU-M 

X1«1.0D0/(2.0D0*(ETS**M) ) 

U( 2, IR)«2 . 0D0*(U( 1, IR)-U( 1 , IR+2 ) )/3 . 0D0 

SK-  U( 2 , IR) *Xl+PH0+PHl 

W<  3, IR)-SK 

LKMM-LKMM-2 

DO  501  K**KFF ,  LKMM ,  KII 

KX-IR+K 

Gl-<K+l)*(K+2.0D0)/( ( 2*K+1 ) * ( 2*K+3 ) ) 
G2-(K+l)*(K+1.0D0)/( ( 2*K+1 ) * ( 2*K+3 ) ) 

G2«G2+K*K/( (2*K-1.0D0)*{2*K+1) ) 

G3«(K-l)*K/( (2*K-1.0D0)*(2*K+1) ) 

U ( 2 , KX ) «U  U , KX ) * ( 1-G2 ) -G1*U< 1 , KX+2 ) -G3*U( 1 , KX-2 ) 

TK«U(2,KX)*(2*K+l)*PE<KX)*Xl 

SK-SK+TK 

W( 3,KX)-SK 

A2-DABS ( TK/SK ) 

IF(A2.LT.ERK.AND.ALT.LT.ERK)  IEE-IEE+1 
ALT-A2 

IF( A2 .GT.A2T.OR. K. LT. LK/4 )  GO  TO  542 
NKT-K 

I F ( A2 . GT . 1 . OE-1 5 )  A2T-A2 
542  CONTINUE 

C*** 

C  DIAGNOSTIC  PRINT  3 
C  WRITE( 21 , 988 )  M,K,TK,SR,A2 
988  FORMAT ( 2I4,3(1PE12.5)) 
r*** 

NKX-K 

IF( IEE.EQ. IEK)  GO  TO  541 

501  CONTINUE 

U(l, IR)«U(2,IR) 

DO  502  K«KFF,LKMM,KII 
KX-IR+K 

U( 1 ,KX)-U( 2,KX) 

502  CONTINUE 
SR-W( 3 , IR+NKT) 

IF( A2T.GT.A2ML)  GO  TO  500 
SKL-SK 
A2ML-A2T 
NKTL-NKT 
500  CONTINUE 
I FLG1--1 
SK-SKL 
NKT-NKTL 
A2T-A2ML 
541  CONTINUE 
RETURN 
END 
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SUBROUTINE  SRUTH ( FTC , LMR ) 

C  SCREENED  RUTHERFORD  LEGENDRE  SCATTERING  COEFFICIENTS 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C*  *  *  * 

COMMON/Cl/IR , IS , JTI , LL , ITRN , LLPl , LLP 2 , ITC , I  PM , LS , LYY ,  JT , JX , 

*  IDFF , LN , LNP , LK , LKP , IEN , I EK , ITR ,  IP , IXY , IPL,ITPE,NS, 

*  KFF , KI I , LM , I  ST , IN , INF , NKT , I FLGl , NKX , IMXM, LAP , LOP , MLT 
COMMON/C2/IDEN( 18 ) , Q0 , AMO , ETA, ALA , BE , Bl , SIO , ERN, ERK , DC , SO , DS , X0 , DX 
COMMON/C  3/P I , X , S , Q , QY , TX , SXl , ET , TS , PHI , FD , DPH , WLPl 

COMMON/C  4 / 1 OT 1 , IOT2 , IOT3 , IOT4 , IOT06 ,MX( 2 ) ,MY( 2 ) , EXMX 
DIMENSION  FTC (LMR) 

FTC ( IR)-1 .0D0 
IF(LL.EQ.O)  RETURN 

FTC ( IR+1 ) -1+ETA-ETA* ( 1 . 0D0+ETA/2 . 0D0 ) *DLOG( 1 . 0D0+2 . ODO/ETA) 
IF(LL.EQ.l)  RETURN 
AX-FTC( IR+1 ) 

FC-DEXP(-EXMX) 

RT-1 . 0D0+ETA-DSQRT( ETA* ( ETA+2 . 0D0 ) ) 

ALP-0 . 25D0* ( 5 . 0D0+RT*RT )/( 1 . 0D0-RT*RT ) 

L1-MLT+LLP2 

ALl-Ll 

WO-FC*RT*DSQRT ( ALl/ALP+1 . 0D0 ) 

W1«FC*DSQRT( ( AL1-1 . 0D0 )/ALP+l . 0D0 ) 

LMX— EXMX/DLOG  (  RT ) 

LSG-LL 

I F ( LMX . GT . Ll )  GO  TO  144 
Ll-LMX 

LSG-L1-MLT-2 
WRITE ( IOTl ,  88  )  LSG 

88  FORMAT ( 'LEGENDRE  SERIES  TRUNCATED  AT  LL-',I4) 

144  CONTINUE 

DO  50  LP-1 , Ll 
L-Ll-LP+1 

W2- ( ( 2*L+3 ) * ( 1+ETA ) *Wl- ( L+l ) *W0 )/( L+2 ) 

FTC ( IR+L)-W2 

WO-Wl 

W1-W2 

50  CONTINUE 
GX-AX/FTC( IR+1 ) 

FTC ( IR+1)-AX 
DO  51  L-2 , LSG 

FTC( IR+L ) — GX*FTC( IR+L) 

51  CONTINUE 
RETURN 
END 

SUBROUTINE  TRNSC( FTC, LLR ) 

:  TRANSPORT  CORRECTED  LEGENDRE 

SCATTERING  COEFFICIENTS 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/Cl/IR, IS, JTI, LL, ITRN, LLPl, LLP2 , ITC, IPM, LS , LXX, JT, JX, 

*  I DFF , LN , LNP , LK , LKP , I EN , I EK , ITR , I P , IXY , I PL , ITPE , NS , 

*  KFF , KI I , LM , 1ST , IN , INF , NKT , I FLGl , NKX , IMXM , LAP , LOP , MLT 
COMMON/C2/IDEN ( 18 ) , QO , AMO , ETA , ALA, BE,Bl,SlO, ERN , ERK , DC , SO , DS , XO , DX 
COMMON/C  3/P I , X , S , Q , QY , TX , SXl , ET , TS , PHI , FD , DPH , WLPl 
COMMON/C4/IOT1 , IOT2 , IOT3 , IOT4 , IOT06 , MX ( 2 ) , MY ( 2 ) , EXMX 

COMMON/C 5/  PHO , ETS , ETl , A2T 
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DIMENSION  FTC(LLR) 

WLPl *FTC ( LLR ) 

Tl-1 . 0D0-WLP1 
IF(LL.LE.l)  GO  TO  11 
DO  10  L-1,LL 

FTC ( IR+L ) ■ ( FTC ( IR+L ) -WLPl )/Tl 

10  CONTINUE 

11  CONTINUE 
ALA-ALA/T1 
LL-LL-1 
LLP1-LL+1 
LLP2-LL+2 
LLR-LLR-1 
RETURN 
END 

SUBROUTINE  DIFF 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

c*  *  *  * 

C  DIFFUSION  CALCULATION 
C  DIFFUSION  COEFFICIENT-DC 

C  ★  *  *  * 

COMMON/Cl/IR , IS , JTI , LL, ITRN , LLPl , LLP2 , I TC , I PM , LS , LX , JT , JX , 

*  IDFF , LN, LNP , LK , LKP , IEN, IEK, ITR, IP, IXY, IPL,ITPE,NS, 

*  KFF, KII , LM, 1ST, IN, INF , NKT , I FLGl , NKX, IMXM , LAP , LOP , MLT 
COMMON/C2/IDEN ( 18 ) , Q0 , AMO , ETA , ALA, BE , Bl , SIO , ERN , ERK , DC , SO , DS , X0 , DX 
COMMON/C3/PI , X , S , Q , QY , TX , SX1 , ET , TS , PHI , FD , DPH , WLPl 

COMMON/C 5/  PHO , ETS , ETl , A2T 
IF(JTI.GT.O)  GO  TO  10 
T1-X*X/(4.0D0*DC*S) 

T2«DEXP(-T1 ) 

T3-4 .0D0*PI*DC*S 

T3-T3**0.5 

DPH-Q0*T2/T3 

IF(IS.EQ.l)  RETURN 

DPH-DPH* ( 1+3 . 0D0*AM0*X/( 2 . 0D0*S ) ) 

RETURN 
10  CONTINUE 

Xl»X/DSQRT( 4 . 0D0*DC*S ) 

CALL  ERF ( XI , ER ) 

DPH-Q0* ( 1 . 0D0-ER )/2 . 0D0 
IF(IS.EQ.l)  RETURN 
T0-3.0D0*AM0*DSQRT(DC/(PI*S) ) 

DPH-DPH+Q0*T0*DEXP(-X1*X1 )/2 . 0D0 

RETURN 

END 

SUBROUTINE  LEG( P , AX , LR , LKR ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

**** 

:  LEGENDRE  POLYNOMIAL  ROUTINE 
:  DETERMINES  LEGENDRE  POLYNOMIALS 
:  FOR  L-0,LK 

l*  *  *  * 

COMMON/Cl/IR,IS, JTI ,LL, ITRN, LLPl, LLP2 , ITC , IPM, LS , LX, JT, JX, 

*  IDFF , LN , LNP , LKX , LKP , I EN , IEK , ITR , I P , IXY , IPL , ITPE , NS , 

*  KFF , KI I , LM , 1ST , IN , INF , NKT , I FLGl , NKX , IMXM , LAP , LOP , MLT 
DIMENSION  P(LKR) 
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P( IR)-1 .0D0 
P(  IR+D-AX 
IF(LK.LE.l)  RETURN 
DO  10  J-2,LK 
JR-IR+J 

PXX-( (2*J-1.0D0) *AX*P( JR-1)-(J-1.0D0)*P( JR- 2 ) )/J 
P( JR)-PXX 
10  CONTINUE 
RETURN 
END 

SUBROUTINE  ERF ( X , ERFF ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  A( 5 ) 

DATA  ITE/-1/ 

IF(ITE.GT.O)  GO  TO  10 
A(l)-0. 254829592 
A(2)— 0.284496736 
A(3)-l. 421413741 
A(4)— 1.453152027 
A( 5)-l. 061405429 
P-0.3275911 
ITE-1 

10  CONTINUE 

T-l . 0D0/( 1 .0D0+P*X) 

ERFF-0 
DO  1  M-1,5 

ERFF-ERFF+A( M ) * ( T*  *M ) 

1  CONTINUE 

ERFF-1 . 0D0-ERFF*DEXP ( -X*X ) 

RETURN 

END 

SUBROUTINE  OUT( AB ,OR , FTC , PHI , AMERR , LA , LO , LLR ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

7  OUTPUT  EDIT  ROUTINE 

COMMON/C1/IR ,  I S ,  JTI  f  LL ,  I TRN ,  LLPl ,  LLP2  f  I TC ,  I  PM ,  LS ,  LX ,  JT ,  JX , 

*  IDFF , LN , LNP ,  LK , LRP , IEN ,  IEK , ITR, IP , IXY ,  IPL , ITPE , NS , 

*  KFF , KI I , LM , I ST , IN , INF , NKT , I FLGl , NKX , IMXM , LAP , LOP , MLT 
COMMON/C2/I DEN (18 ) , Q0 , AMO , ETA , ALA , BE , Bl , SI 0 , ERN , ERK , DC , SO , DS , X0 , DX 
COMMON/C 3 /P I ,  X ,  S ,  Q ,  QY ,  TX ,  SXl ,  ET ,  TS  ,  PHI ,  FD ,  DPH , '»7T.?1 

COMMON/C 4/IOTl , IOT2 , IOT3 , IOT4 , IOT06 , MX( 2 ) ,MY( 2 ) , EXMX 
DIMENSION  PHI ( LA, LO) , AMERR( LA, LO ) , AB ( LA ) ,OR(LO) 

*  , FTC ( LLR ) 

REAL  FS( 3 ) , FL( 5 ) , YC,XC 
DATA  ITU/1/ 

WRITE( IOTl , 399 )  IMXM 

399  FORMAT( //IX, 'MINIMUM  CONTAINER  ARRAY  DIMENSION- ', 16 ) 

IF(IPL.NE.l)  GO  TO  400 
PHMX-0.0D0 
DO  401  1-1 , LAP 
DO  402  J-l , LOP 

IF(PHI( I , J) .EQ.7 .7777E17 )  GO  TO  402 
I F ( DABS ( PHI ( I , J ) ) . GT . PHMX )  PHMX-DABS ( PHI ( I , J ) ) 

402  CONTINUE 
401  CONTINUE 
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400  CONTINUE 

IF(JTI.EQ.l)  GO  TO  16 
FL(1)-4HAT  X 
FL(2)-2H-0 
IJ-2 

GO  TO  14 

16  CONTINUE 

FL ( 1 )-4HCONS 
FL(2)«4HTANT 
FL ( 3  )«4H  HAL 
FL ( 4 )  -4HF-SP 
FL( 5)-3HACE 
IJ-5 

14  CONTINUE 

I F ( IS.EQ.2 )  GO  TO  18 
FS( 1 )«4HISOT 
FS ( 2 ) -4HROPI 
FS( 3)«1HC 
GO  TO  15 
18  CONTINUE 

FS( 1 )-4HANlS 
FS ( 2 ) —  4  HOTRO 
FS(3)-3HPIC 

15  CONTINUE 

I F ( MX ( 1 ) . EQ . 1 )  XC-lHS 
I F ( MX ( 2 ) . EQ . 1 )  XC-1HX 
I F ( MY ( 1 ) . EQ . 1 )  YC-4H  S/ 

IF(MY<2) .EQ.l)  YC-4H  X/ 

17  CONTINUE 
NT- LAP/6 
LTX-5 
IUX-6 

IF(ITR.GT.O)  GO  TO  36 
NT-LAP/1 1 
LTX-10 
IUX-11 
36  CONTINUE 

NX-LAP-IUX*NT 
IF(NX.GT.O)  NT-NT+1 
IF(ITU.LT.Q)  GO  TO  333 
WRITE( IOT1 ,222) 

222  FORMAT( /IX, 'PLANE  GEOMETRY  '//IX, ' SOURCE :' ) 
WRITE( IOTl , 23 )  ( FL{ I ) , 1-1 , I J ) 

23  FORMAT ( 3X , 6A ) 

WRITE( IOTl, 230)  (FS (I), 1-1,3) 

230  FORMAT ( 3X , 6A ) 

IF( IS. EQ. 2 )  WRITE( IOTl , 38 )  AMO 
38  FORMAT ( 3X, 'MU0-' ,1PE12. 4) 

WRITE( IOTl, 66)  LL, ( FTC ( LY ) , LY-IR, LLR ) 

66  FORMAT( /IX, 'SCATTERING  ORDER  LL-',I4/ 

*  '  LEGENDRE  COEFFICIENTS:'/ 

*  ( IX, 5( 1PE12 . 5 ) ) ) 

IF(ITC.GT.O)  GO  TO  84 
WRITE ( IOTl , 85 )  ETA, ALA 

85  FORMAT ( /IX , ' ETA- ' ,1PE12.4,2X, 'ALA-' ,1PE12.4) 
GO  TO  86 
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84  CONTINUE 

WRITE ( IOTl , 87 )  ALA 
87  FORMAT ( /IX , 'ALA- ' , 1PE1 2 . 4 ) 

86  CONTINUE 

WRITE( IOTl , 81 )  BE,SI0 

81  FORMAT ( /IX , '  BE  - ' , 1PE1 2 . 5 , 2X , ' SI 0- ' , 1PE12 . 5 ) 

IF< IDFF.LT.O)  WRITE( IOTl ,  65  ) 

65  FORMAT( /IX, 'DIFFUSION  THEORY  COMPARISON') 

IF( IDFF.GT.O)  WRITE ( IOTl , 74 ) 

74  F0RMAT(/1X, ' FLUX  CALCULATION') 

333  CONTINUE 

DO  1  LLT“1 , NT 
I F ( ITR.LT.O) 

+WRI TE ( IOTl , 32 ) 

32  FORMAT ( / ) 

IF(XPL.EQ.l)  WRITE ( IOT2 , 34 )  LLT, LAP , LOP 

34  F0RMAT(1X, 'FILE' ,12,214) 

IXl“l+IUX* ( LLT-1 ) 

IF(NX.EQ.O)  GO  TO  35 
IF(LLT.EQ.NT)  LTX-NX-1 

35  CONTINUE 
IX2-LTX+IX1 

WRITE( IOTl , 20 )  YC , XC ,  < AB ( IX ) , IX-IX1 , 1X2 ) 

IF( IPL.EQ.l ) 

*WRITE ( IOT2 , 200 )  YC,XC,  < AB ( IX ) , IX-IXl , 1X2 ) 

C*  *  *  * 

C  FORMAT  MAY  BE  CHANGED  TO  ACCOMMODATE  PLOTTING  PACKAGE 
200  FORMAT ( 2X, 2A4 , IX, 6 ( 1PE11 . 4 ) ) 

"*  *  *  * 

20  FORMAT! /2X, A, A, 1X,10(1PE11.4) ) 

I F ( ITR.GT.0) 

♦WRITE! IOTl , 21 ) 

21  FORMAT ( IX, 38 ( 2H**  ) ,  1H* ) 

IF(ITR.LT.O) 

♦WRITE! IOTl, 22) 

22  FORMAT( IX, 61 ( 2H** ) ) 

DO  30  IU-l,LOP 

WRITE ( IOTl , 31 )  OR( IU ) , ( PHI ( IY, IU) , IY-IXl , 1X2 ) 

31  FORMAT( 1PE11 . 4 , 11  ( 1PE11 . 4 ) ) 

IF(IPL.NE.l)  GO  TO  30 
DO  300  ICK»IX1,IX2 

I F ( PHI ( ICK , IU ) . EQ .7.7777E17)  PHI ( ICK, IU)-5 . 0D0*PHMX 
300  CONTINUE 

WRITE! IOT2, 331)  OR(IU) , ( PHI ( I Y, IU) , I Y-IXl , 1X2 ) 

:  FORMAT  MAY  BE  CHANGED  TO  ACCOMMODATE  PLOTTING  PACKAGE 
331  FORMAT! 1PE11 . 4 , 6 ( lPEll . 4 ) ) 

*  *  *  * 

30  CONTINUE 

IF! ITR.GT.0) 

+WRITE(IOTl,21) 

IF( ITR.LT.O) 

+WRITE(I0T1,22) 

WRITE! IOTl, 37) 

37  FORMAT! //6X, 'NUMBER  OF  TERMS  FOR  CONVERGENCE') 

DO  40  IU-l,LOP 
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WRITE ( IOTl ,33)  { AMERR( I Y , IU ) ,IY-IX1,IX2) 
33  FORMAT( 6X, 11 ( 2PE12 . 5 ) ) 

40  CONTINUE 
1  CONTINUE 
ITU— 1 
RETURN 
END 
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Abstract 


This  document  describes  the  operation  of  the  MSLAB  program  which  embodies  the 
F/v  algorithm  for  1-D  neutral  particle  transport  in  heterogeneous  slab  geometry.  The 
F/v  algorithm  was  first  developed  by  C.  E.  Siewert  and  has  been  shown  to  be  one 
of  the  most  robust  of  any  of  the  1-D  transport  methods.  The  theory  and  numerical 
development  of  the  algorithm  as  well  as  programming  notes  and  sample  problems  are 
included. 


I  Introduction 


The  determination  of  the  dose  distribution  resulting  from  electron  energy  deposition  along 
the  path  of  travel  in  some  applications  requires  the  solution  of  the  neutral  particle  trans¬ 
port  equation.  Currently,  there  are  several  numerical  solution  methodologies  that  provide 
reliable  solutions  to  the  electron  transport  equation  including  Monte  Carlo,  eigenfunction 
expansions,  P/v  theory,  numerical  transform  inversion,  and  S/v  theory.  The  last  of  these, 
discrete  ordinates  or  S/v  theory,  has  proven  to  be  a  particularly  valuable  approach  since 
it  is  more  cost  effective  than  most  others.  However,  even  with  its  efficiency  and  appar¬ 
ent  accuracy,  the  S/v  method  requires  complete  discretization  of  the  independent  variables 
(space  and  direction)  leading  to  the  question  of  just  how  reliable  is  the  S/v  solution.  For 
this  reason,  it  is  imperative  that  highly  accurate  independant  solutions  continue  to  be 
developed  in  order  to  assess  the  numerical  discretization  error  inherent  in  the  discrete 
ordinate  approximation. 

'Work  performed  for  the  CUBUC  contract  (No  S-0-7542)  for  Rome  Air  Development  Center,  Hanscom 
Air  Force  Base. 


This  manual  gives  the  details  of  one  such  numerical  monoenergetic,  1-D  transport  deter¬ 
mination,  the  F/v  method,  which  was  originally  developed  by  C.  E.  Siewert  [1]  but  was 
never  released  as  an  operational  code.  The  major  advantage  of  the  F jv  method  is  that  only 
a  single  approximation  (in  the  angular  variable)  need  be  made  allowing  the  independant 
variables  to  be  treated  as  continuous. 

This  manual  will  delineate  in  detail  the  mathematical  derivation,  numerical  analysis,  pro¬ 
gramming,  and  the  usage  of  the  resultant  program.  Section  II  provides  the  mathematical 
analysis  required  to  generate  the  coupled  equations  from  the  transport  equation  that  lead 
to  the  boundary  fluxes.  In  Section  III  the  F/v  approximation  for  a  single  slab  and  multiple 
slabs  with  isotropic  sources  is  given.  Section  IV  details  the  numerical  analysis  required  to 
implement  the  F/v  method  in  FORTRAN  77.  The  last  section  details  the  operation  of  the 
code  MSLAB  including  a  flow  chart,  explanation  of  the  input  deck  plus  several  examples, 
explanation  of  the  output  generated,  and  output  from  the  examples.  In  addition,  detailed 
programming  notes  deemed  necessary  by  the  author  are  provided. 


II  Coupled  Integral  Equations  for  the  Exiting  Flux 


With  only  minor  differences,  the  theory  associated  with  the  F.v  method  applied  here,  is 
contained  in  [1],  For  clarity,  the  major  points  of  the  theory  will  be  summarized. 

The  one-group  neutral  particle  transport  equation  to  be  solved  in  a  slab  of  width  A  for 
|// 1  <  1  and  a  <  x  <  b  is 


d  i 


^  $(1-,//)  -I-  ^S(r) 


(la) 


where  isotropic  scattering  and  an  isotropically  emitting  spatial  source  have  been  assumed. 
The  boundary  conditions  are 


$(«,/*)  =  FL(/i),  //  >  0 


(lb) 


$(fc./z)  =  F n(-/i).  //  <  0 


(lc) 


with  F|,  and  F|<  general  functions  of  //.  These  functions  are  assumed  known  but  they  may 
actually  depend  upon  the  transmitted  and  reflected  fluxes  from  adjacent  slabs  as  indicated 
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below.  Indeed,  it  will  be  the  boundary  fluxes  that  connect  all  slabs  together  as  shown  in 
Section  III.B. 


As  a  starting  point,  the  i’,v  method  requires  integral  equations  for  the  reflected  and  trans¬ 
mitted  fluxes  from  the  slab.  These  equations  are  derived  by  manipulation  of  eqs.  (1)  as 
follows: 


1.  Let  fi  be  replaced  by  — p,  multiply  by  e  T/s  and  integrate  over  r  on  [aj>]  to  give 


-^-B  in,s)~  /dr  $(r,  -n)e  »  =  —  />*(«)  +  —  S*(s,  a,  b) 

H  —  S  Ja  l  fl  —  S  Z  fl  —  S 


(2a) 


where, 


_£ 


B(//,s)  =  $(a,  -fi)e  > 


(2b) 


>*(-■0  =  J‘-~e  *  j 


(2c) 


S"(.s,«,6)=  [dr  S(r)e~  * 

J  a 


(2d) 


and  s  is  a  complex  variable  in  the  cut  plane  [—1,1]. 
2.  Integrate  eq.  (2a)  over  fi  on  [—1, 1]  to  obtain 


y"+1,  /*  \  M*)pm{s) 

/  d// - B(/i,s)  = - 

J- 1  fl  —  $  s 


+  L(s)S‘(s,  o,  6) 


(3a) 


where, 


(3b) 


e.s  /•+'  1  cs  ,  /  s  —  1 


(3.-: 


3.  Multiply  eq.  (3a)  by  and  restrict  .s  to  Re(.s)  >  0  to  give 


/. 


d//  — - — C’( // ,  s)  =  A(s) — — -  -f  L(s)S’(.s,  «. 

i  //  —  . S  .5 


(4a) 
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where 


fb  /+1 

T(s)  =  J  dr  e  »  j  d,  $(r,  —  fx) 


C(fi,s)  =  $(a,  ~n)  -  $(6,  ~n)e  * . 


(4b) 

(4c) 


4.  To  obtain  a  second  equation,  replace  s  by  —  s  and  fx  by  —fx  in  eq.  (3a)  and  multiply 
it  by  e-fc/s,  again  with  Re(s)  >  0  to  give 


/  dfi — — — D(/i,  s)  =  A(s) — — -  -f  L(s)S*(  —  s,  a,  b) 
J-i  LI  —  s  s 


6 

)e  '* 


where 


r  b-r  r+l 

J’(s)  =  J  dr  e  *  j  dfi  $(r,  ~fx) 

D(//.,  s)  =  $(/>,//)  -  $(a,/i)e~^. 

5.  In  addition,  since 

A(±i/0)  =  0,  v0  >  1 

when  s  —  v0  in  eq.  (4a)  and  eq.  (5a),  we  have 
r+i 


CV 0  /,+  ‘ 


i  — - - C(/x,  r'o)  =  ~S*(ir0,  a,  &)e‘ 

H  -  v0  2 


[  d/i  — — — D(/i,//0)  =  ~S*(-i/0,«,6) 

2  J-i  //.  —  v0  l 


_ fc. 

t  vo 


(5a) 


(5b) 

(5c) 

(6) 


(7a) 

(7b) 


Equations  (4a)  and  (5a)  are  coupled  integral  equations  with  constraints  given  by  eqs.  (7). 
These  integral  equations  will  now  be  restricted  to  s  on  the  cut  [0, 1].  This  is  accomplished 
using  the  Plemelj  relations  in  the  form 


n  in 


(  1  )-r 

1 

V  //  -  (v  ±  if )/ 

±  —  v } 


where  //,  f.  [().  1]  and  P  [  ]  represents  a  principal  value  when  under  an  integral.  Thus,  if 
\(s)  is  defined  as  the  integral 
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/+» 

&n 


+\  fit*) 


H  —  s 


then 


Inn  [xiu  ±  ie)\  s  xV)  =  vj [dp  ±  i*f(v). 


(8) 


After  application  of  eq.  (8)  to  eq.  (4a),  we  find 


V  f+ld,t  _JL_C(p,  v)  i  t7ri/C(i/,  t/)  =  A  a,  b)e* 

J-\  p-v  v 


(9a) 


noting  that  I*  and  S*  are  analytic  functions  and 


cv 


A±(i/)  =  A(n)  ±  <7r  — 


(9b) 


cv  . 

A(i/)  =  1  +  —  In 


1  +  v 


1  -  v 


(9c) 


m 


L±(i/)  =  L{v)  ±  — 


(9d) 


L  ( iv )  =  —  In 


1  -  iv 


1  +  v 


(9e) 


Next,  the  plus  and  minus  branches  from  eq.  (9a)  are  subtracted  to  ghe 


u 

c  " 


=  -C(v,v)  -  — S*(i '.«,&); 
v  r  cv 


(10a) 


the  plus  and  minus  branches  of  eq.  (9a)  axe  then  added  to  give 


■p[+du  —¥—C (/vv,  =  +  L(  I/)S*(  v,a./> )r*. 

J- 1  fi  ~  v  i' 


(10b) 
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Substitution  of  eq.  ( 10a)  into  eq.  ( 10b)  then  yields  the  first  of  the  desired  integral  equations 


cis  r+i  u  p'Z 

u\{i/)C{u,u)  —  —  V  dp - C(p,z/)  +  —  S*(i/,  a,b).  (11a) 

1  J- i  p  —  v  2 

Similarly,  for  eq.  (5a),  we  obtain  the  second  integral  equation 

u\(v)D{v,is)  =  ^  vf  dp  —  -D(p,t/)  +  €—  S*(-i/,a,6).  (lib) 

2  7-i  p  —  v  2 

Equations  (11)  are  the  integral  equations  to  be  solved  subject  to  the  constraints  given  by 
eqs.  (7).  The  Fat  approximation  is  now  constructed. 

Ill  The  F,y  Approximation 

III. A  Single  Slab 

Equations  (11)  can  be  written  in  terms  of  the  reflected  and  transmitted  fluxes  as  follows: 
A (//)$(«,-»')  ~  (~vf  d/i  — - — $(a,~p)+ 

2  io  p  —  v 

+  /  d/i  -$(fr,/i)  =  ~  /dp  — FL(/i)  +  (12a) 

2  Jo  fi  +  v  l  Jo  p  +  v 

+(  (  A(t/)Fn(ii)  -  ^  Vf  d/i  Fr(/z)1  -1-  ^-S*(i/,  a,b) 

(  2  7o  p  —  v  )  2u 

\{is)${b,i>)  -  -  P / d/i  — d>(A,  p)-f 

2  7»  /i  —  j/ 

+  /<l/'  ~ — $(«• -/')  =  -:  /d/i  — 7— Fu(/i)  +  (12b) 

fi  |  A(  ii)F|,(  fi)  -  ^  vf  d/i  — - — F|,(/i)l  4-  —  S*(  —i/,  a,  A). 

(  2  A>  /i  —  //  J  2/^ 

The  corresponding  equations  for  the  constraints,  eqs.  (7),  are  obtained  from  equations  (12) 
by  letting  u  be  and  A(i/0)  be  zero. 
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The  F/v  approximation  is  simply 


$(a,  -n)  =  FR(p)e  *  4-  ^  £  a*^o,(n) 

“  a=0 

(13a) 

=  FL(A*)e~^  +  1  S 

^  n 

(13b) 

*  a=0 


where  aa,  and  ba  are  coupling  coefficients  which  are  determined  via  collocation  as  will 
be  shown;  and  0  are  a  set  of  known  basis  functions.  When  the  F/v  approximation  is 
introduced  into  eqs.  (12),  we  obtain 


N-l 

[«0B0(i/)  4-  f>0Aa(i/)e_^j  =  Rj(i/,a,6,  A) 


a=0 


(14a) 


N-i 

J2  [floA0(i/)e’^  +  6QBQ(r')j  =  R2(i/,a,6,  A) 

Or =0 


(14b) 


where 


Aa(i')  =  £  f  d/l  -  2 — I /’,-,(/<) 

2  Jo  n  4-  V 

Bl,('/)  =  A(  »>-)»/’.,(")  -  7>  P\ d/‘  — t/’oOO 

2  Ju  fi  —  v 


(14c) 


( 14d) 


Ri(^,«,  b,  A)  =  Jdfi  ft  [S(  A,  //,  v)F iXji)  4-  C(A,p,  i/)Fr(//)]  4- 


(14c) 


+  —  Sm(v,a,b) 
cv 


R2(i/,«,6,  A)  =  Jd/i  fi  (S(A,/x,i/)F»(p)  4-  C(A,/x,i/)FL(/z)j  4- 


(14f) 


c  » 


4 - S'(  — //,  a,  b) 

cv 
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C(A ,/i,  u) 


-A 

e  <*  —  e 


[X  —  v 


(14g) 


S(A,/i,i/) 


/x  +  1/ 


(14h) 


The  coefficients  aa  and  ba  are  obtained  using  collocation  where  1/  is  specified  to  be  j/0  plus 
the  N  —  1  zeros  of  the  shifted  Legendre  polynomials,  Pn_i(i/),  and  the  roots  being  up, 

0  =  1,2,3, ...  ,N. 

Note  that  for  u  =  i^o,  eqs.  (14)  become  the  constraints  specified  by  eqs.  (7)  [with  the 
convention  that  A(i/0)  =  0).  The  2N  linear  equations  to  be  solved  for  aa  and  ba,  a  = 
0, 1, 2, .  . . ,  N  —  1  are  therefore 


N  —  I 

E 

i»=0 


«J.,a  +  l>„Aaiic 


RUi(a,b,A) 


(15a) 


N  - 1 

E 

<.=o 


a0A(,pt  **  +baBap 


R2p{a,b,A) 


(15b) 


where 


(15c) 

B,»^  —  Bft(  ) 

( 1 5d ) 

r3 

i> 

ii 

» 

> 

(15c) 

R  2p(a,h,A)  =  R2(u.i.(i.b.A). 

(15f) 

From  experience,  it  has  been  found  that.  eqs.  (13),  for  the  reflected  and  transmitted  fluxes 
respectively,  do  not.  provide  exceptionally  accurate  results,  especially  near  //  =  0.  Thus 
C.  E.  Siewert  in  [1]  has  developed  a  post  processor  which  provides  a  more  accurate  evalu¬ 
ation.  This  expression  is  derived  with  minor  modification  by  noting  that  from  eq.  (14d) 
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BaM  =  A*(l/)j/>a(l>)  -  Bma(u) 


(16a) 


where 


A’(i')  =  1  +  —  In 


is 

1  + is 


(16b) 


B-(iz)  =  r 


2 


llj'aiH)  -  VlJ>0(v)' 

(l-  is 


(16c) 


When  eq.  (16a)  is  introduced  into  eqs.  (14),  there  results 


4>(a, -/i)  =  Fn(/t)c  >-  + 


(17a) 


N- 1 


+  5~A~(/7)  +  Y1  [aaB‘(/<)  -  baAmJfi,)} 


d>(/),/0  =  Fl(/0<-~  “  + 

+  0X'(ft)  |R2 (/q«A  A)  +  (6->B;(/x)  -  rt(tAt;(/i)j 


(17b) 


with 


A*  (// )  =  A,, (//  )t  '  >• . 


(17c) 


Singular  integral  equations  can  also  be  obtained  for  the  flux  at  a  position  .»•  in  the  slab 
interior  by  first  letting  a  =  .r  in  eq.  (12a)  to  give 


(18a) 


BIO 


and  letting  b  —  x  in  eq.  (12b)  to  give 


A(j/)$(x,i/)  -  |  vld^x  — —  ${x,fi)  + 

2  Jo  [x  —  v 

C  j-a  ix  c  ['  IX 

+-e  -  dfx——-$(a,-ix)  =  -  dn — - — $(x,-/x)  + 

2  Jo  ix  +  v  2  Jo  u  +  v 


+e-V  {AW)FL(^)  FL(,i)}  +  ^S-(-^,a,x). 


Then,  if  the  approximations 


$(x,  -/0  =  F R(fx)e  “  +  ^  51 


N-i 


<>=0 


N-I 


$(.!■,//)  =  F|,(/l)t  '•  +  -  5Z  4><M/0 


a=0 


are  substituted  into  eqs.  (18)  and  the  result  evaluated  at  e  =  up,  the  following  line 
of  equations  to  be  solved  for  ca  and  da  are  obtained 


N-l  „  _t_xN-l 

(c«B«/3  “  rfoAa/(]  =  Ri^(a,  b ,  x)  -  c  v»  J2  baAap 

o=0  o=0 


N-l 

5  ^  [  ^*o "4"  —  ri'2/d A,  x ) 

o=0 


when* 


_£^a  N~' 

e  ^  5Z  a«Af.;) 

o=0 


Ri(  //,  a,  /).  ./• ) 


/'  [Fi.(/»)S(fc  -  x,n,v)<  <■  +  Fn(/i)C(/>  -  .i-,/i,t/)]  + 

_r 

4 - S"(".  .r.  />) 

ru 
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(18b) 


(19a) 


(19b) 

ar  set 


(20a) 


(20b) 


(20c) 


R-2 (v,a,b,x)  =  Fl(/x)C(x  -  a,  n,  u)  +  Fr(h)S(x  -  a,  /x,  u)e~  »  +  (20d) 


_  -JL 

+- — S*(  —v,  a,  x ) 
cv 


Rj ,2ffW,  b,  x)  =  Rli2(j//),a,  b,x). 


(20e) 


Since  aa  and  ba  have  previously  been  determined,  this  set  of  equations  can  be  solved  for 
ca  and  da.  The  post  processing  step  is 


$(.r,  -p)  =  Fu(/0e  <*  + 


+  TV^I)  +  Y1  {coBa(d)  +  (da  -  bae  m  'j  Aa(/x)| 


(21a) 


$(.*■,//)  =  FL(/i)r  -  + 


c  r  n-i 

+  ~ o,fc,;r)  +  52  {rfoB;(/0  +  (c0  -  aQ<T“)  Aa(/i 

00  ,t=o 


(21b) 


III.B  Multiple  Slabs 

So  far,  only  a  single?  slab  has  been  considered  We  now  consider  a  heterogeneous  medium 
composed  of  a  set  of  slabs  bordering  a  vacuum.  On  the  left  vacuuin/slab  interface,  a  beam 
source  of  strength  So  and  direction  fi0  is  assumed  to  impinge.  On  the  right  vacuum/slab 
interface,  there  is  no  source.  The  connection  between  each  slab  is  made  through  the 
boundary  term  on  the  right  hand  side’s  of  the  eqs.  (15)  and  eqs.  (20).  Three  distinct  cases 
can  be  identified  including: 

Lvs-Rss:  left  vacuum  with  a  beam  source/slab  interface,  right  slab/slab  interface. 
Lss-Rss:  left  and  right  both  slab/slab  interface,  and 
Lss-RsV:  left  slab/slab  interface,  light  slab/vacuum  interface. 


The  particular  right  hand  side’s  for  each  interface  will  now  be  specified.  Since  R  i  ( a,  />,  A )  = 
R  i( t-\ a.  b,  u)  and  R  2(i'.«Ji,A)  =  R2(u.  a,  h,  b).  only  R|  and  R2  need  be  considered  in  the 
following  description. 
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1.  Lvs-Rss 

For  this  case,  the  slab  is  the  initial  slab  going  from  left  to  right  and 

Fl(^)  =  S0 6(n  -  fi0),  /A/r o  >  0 

Fr  =  $2(Ai,  —  M)i  ll  >  0 

where  Ai  is  the  thickness  of  the  first  slab  and  $2(Ai,  — //)  is  the  reflected  flux  from 
the  second  slab.  Then  from  eqs.  (20c, d)  (for  a  —  0  and  b  =  Aj) 


Ri(i i/, 0,  Ai ,  x)  —  fi0 S0S( Ai  -  x ,  (J.Q,  v)e~*>  4-  ^S’(i/,  x ’ 

+  dn  fxC( A,  -  x,/i,v)$2(Ai,-/i)  -  e  u» 

Jo  O=o 


(22a) 


R2(r/,  0,  Ai,x)  =  /.i0S0C(x,fi0,is)  4-  -^-S*(  — u,0,  x)+ 

i  a  -i  N-1 

4  / d//  pS(.r,p,  iv)$2(Ai,  -n)e~  -  e  v»  YL  ba kQ0. 

•'O  a=0 


(22b) 


2.  Lss-Rss 

For  slab  k,  the  entering  fluxes  are  (for  a  —  .c jt- i  and  ft  =  Xk) 

_  *k-  1 

Fi.(/i)  =  $*-i(a-*_i,/j)  +  Sof  "°  H o) 

Fr(m)  =  $*+i  (**,-/*) 

where  the  uncollided  component  has  been  explicitly  accounted  for.  Thus 


„  _ x.  C  * 

=  fi0S0S{xk  -  x,no,v)e  +  —Sm{v,x,xk)+ 


4  Jdfi  fi  [S(xfc_,  -  x,/x,i/)$fc_i(xfe_i,^)e  "  4 

_£i-x  N-  1 

+  C(.r k  -  x,ii,v)$k+  \{*k,  ~/i)]  -e  1,0  J2  b°A<*P 

a=0 


(22c) 


ik-i  e  *- 

R2( .rx—  i . )  =  //0S0C(.r  -  xfc-i, ")e  "°  d - S*( -/a xk-\ ,  •?' )4 


^d//p  S(, 


j  jt_i ,  //,  t/)^fc+)(  ,  *,  ~/i)t  "  4 


(22d) 


4  C’(.r  -  j'/t-  \ .  (i .  i*  ~  i  ( xk- 1 ,  /* )]  -  f- 


Y  °<'A< 
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3.  Lss-Rsv 

This  is  the  case  for  the  right  most  slab,  designated  slab  I<,  where 


J  A- 


FL(/i)  =  $A_i(xA_,,p)  +  s0e  6{n-no) 


Fr(//)  =  0 


Therefore  (for  a  =  xA-_!  and  b  =  x A ) 


£ 

e*- 


R,(i/,:rA_i,xA,x)  =  ^0S0S(xk  -  x,p0,  i^)e  H - S*^,  x,xA  )+ 


+  jfdp  pS(xA-  -  X.^l^K-l^/C-l,^) 


cv 

e  **  — 


xA-x  N-  1 

-6  ^  ^  i0Aft|3 

a=0 


J  /\  —  1 


R2(^,  l'K_i,.rK,x)  =  p0SoC(x  —  XA-.i,/i0,  ^  "°  + 

+  - - S*(-i/,XK-i,ar)  +  / dp  pC(x  -  xA  _,,  p,  i/)4>a  _i(xa  _,  , /x) 

ci/  Jo 


J-*/i-i  N  - 1 

“t"  ^  Yla»  Aatf- 

c»=0 


III.C  Uniform  Source 

A  useful  source  distribution  to  assume  is  one  that  varies  exponentially  across 
S(x)  =  S r-e_‘u  r.  Thus,  the  required  S*  for  slab  k  is  simply  given  by 


Sm(v,a,b)  =  S k  f dr  c~(*k+1^T 

J  a 

1  +  m-n  L 

The  appropriate  source  contributions  for  Rj,  R 2.  Rj,  and  R-2  are  for  R2: 


—  S*(  v,  a ,  fr) 


rr 


Sfcr-*" 
c( 1  +  otr) 


[}  _  f  -('U  +  'i)A‘' 


(22e) 


(22f) 


each  slab, 


(23a) 
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for  R2: 


e  * 


cv 


S'{-v,a,b) 


Sbc'“ka 

c(  1  -  CVjfet') 


—  e 


-^k^k 


(23b) 


for  Rj: 


€  y 

— S"(*/,X,Xt) 


cv 


Ske~°k' 
c(l  +  akv) 


1  _  e-(^*-x)(o*  +  i) 


(23c) 


for  R-2- 


€  » 


CV 


S’(v,xk_i 


Ske~akIk~l 
c(l  -  akv) 


*-*k-\ 

e  * 


(23d) 


IV  Numerical  Implementation  and  Demonstration 


IV. A  Numerical  Considerations 


Several  decisions  concerning  the  numerical  method  associated  with  the  F/v  method  must 
be  made,  including: 


•  numerical  quadrature 

•  zero  search 

•  choice  of  vp 

•  choice  of  basis  functions 

•  matrix  inversion 

•  iteration  strategy. 


All  integrals  required  for  the  matrix  dements  and  inhomogeneous  terms  will  be  performed 
using  a  shifted  Legendre-Gauss  quadrature  scheme.  While  the  numerical  integration  is 
straight  forward  for  A,,,  Rt,  and  Rj.  it  is  not  for  B,(.  This  is  a  consequence  of  the  principle 
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value  integration.  However,  upon  regularization  of  the  integral  as  performed  for  the  post 
processing  step,  we  need  only  to  evaluate  B*(iz)  given  by  eq.  (16c)  which  no  longer  presents 
a  difficulty.  In  addition,  to  allow  for  some  flexibility  and  numerical  adjustment,  a  different 
quadrature  order  than  for  the  matrix  elements  is  allowed  for  the  integrals  over  the  boundary 
fluxes  contained  in  the  inhomogeneous  terms  connecting  the  slabs  (Ri  and  R2).  It  should 
be  noted  that  the  evaluation  of  the  matrix  elements  differ  from  the  recurrence  relations 
used  originally  in  [1].  Comparable  accuracy  using  quadratures  however  can  be  achieved. 

A  zero  search  is  performed  to  determine  uQ  from  A(r>o)  =  0.  A  simple  bisection  procedure 
with  zero  bracketing  and  refinement  is  employed. 

The  most  convenient  choice  for  (3  =  1,2,  ...,N  are  the  N  —  1  zeros  of  the  shifted 
Legendre  polynomials  P)(/_1(i')  on  [0, 1].  Other  choices  are  also  available  such  as  uniform 
spacing  and  the  zeros  of  shifted  Chebyshev  polymomials  however  no  significant 

difference  has  been  observed  with  these  choices.  The  choice  of  the  basis  functions  can 
be  crucial  to  the  proper  operation  of  the  algorithm.  I11  the  past,  both  the  monomials 
and  shifted  Legendre  polynomials  P *(/x)  have  been  used.  Because  of  their  ill-conditioned 
behavior,  however,  monomials  have  been  abandoned.  I11  addition  to  the  above  choices, 
several  others  are  possible,  including  shifted  monomials  {2fP  —  l)11  and  modified  Legendre 
polynomials  Pa(2 fP  —  1)  where  7  is  a  free  parameter.  From  experience,  best  results  have 
been  found  for  y  %  0.75. 

The  linear  system  of  equations  generated  by  the  F/v  algorithm  for  both  the  boundary  and 
interior  fluxes  are  generally  of  the  form  (ji  =  0.  1, . . . ,  Ar  —  1 ) 


N  - 1 

7:  [('..A.,/*  +  faHaj]  =  \ i,( 

a=0 


(24a) 


,v- 1 

y*.  +  /„  A„yij  =  \2/i- 

>.=0 


(24b) 


Adding  and  subtracting  these  equations  yields 


,v  -1 


y  )(  a,.,v  +  //,,.) )  j  —  \  ui  +  \2  a 


(25a) 


v- 1 


y  [(f  —  ./.>  )(  A,,,j  —  —  \l.<  —  \  i.1- 


(25b) 
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These  equations  are  solved  for  gf  =  ea  +  fa  and  ga  =  ea  —  fa  using  the  LU  decomposition 
[2],  then 


\  {at  +  da  )  (26a) 

fa  =  \  (gt  ~  9*)  (26b) 


The  following  iteration  strategy'  has  been  adopted  in  order  to  ensure  convergence  of  the 
angular  flux.  The  order  N  of  the  approximation  is  increased  until  the  fluxes  on  the  slab 
boundaries  have  converged  to  a  desired  relative  error.  For  more  than  one  slab,  an  additional 
inner  iteration  at  each  N  is  imposed  since  each  slab  is  considered  separately  and  the 
entering  fluxes  are  not  known  a  priori.  The  inner  iteration,  essentially  additional  cycles 
through  the  slabs,  are  preformed  with  the  boundary  fluxes  updated  after  each  cycle.  TToon 
convergence  of  the  boundary  fluxes,  the  (N)  iteration  for  the  interior  fluxes  is  begun  at 
the  order  of  convergence  of  the  boundary  fluxes  in  order  to  avoid  unproductive  lower  (N) 
order  iterations. 

In  addition  to  the  angular  flux,  the  reflection  and  transmission  coefficients, 

A*  =  ^dp  fi(t>(a, -p)/S0  (27a) 

B*  =  J dp i  ft<t>(b,n)/S0,  (27b) 


are  desired.  Since  the  angular  flux  at  the  boundaries  (a  and  b)  are  known  at  Gauss- Legendre 
quadrature  points,  the  integrations  in  eqs.  (27)  can  be  easily  performed. 


IV.  B  Demonstration 


In  order  to  test  the  accuracy  of  the  M  S 1. A H  code  and  test  the  algorithms  for  proper  coding, 
several  trials  have  been  performed,  the  results  of  which  will  now  be  reported. 

Trial  1  involves  the  verity  at  ion  of  the  accuracy  of  a  single  slab  of  width  1  and  c.  =  0.9  by 
setting  the  quadrature  order  to  30  and  varying  the  desired  relative  error  e.  Table  1  shows 
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the  reflected  and  transmitted  fluxes  for  e  =  10~2,  10-3,  10-4,  10~5  with  the  results  for 
the  last  case  identical  to  those  of  [lj.  As  can  be  observed,  the  results  are  always  within 
the  prescribed  accuracy  giving  confidence  in  the  coding  of  the  algorithm.  As  e  is  further 
reduced  to  10~8,  the  results  remain  invariant.  An  Fyv  approximation  of  N  =  22  was  required 
to  achieve  an  accuracy  of  10-5.  Table  2  is  similar  to  Table  1  for  e  =  10-2,  10-4,  10-6  for 
the  intern  r  point  x  —  0.5.  At  e  —  10-6  the  results  are  identical  to  those  of  [1], 

In  the  second  trial  the  quadrature  order,  Lm,  was  varied  for  the  same  problem  with  e  = 
10~5.  For  Lm  =  16,  the  fluxes  have  converged  to  the  results  of  [1]  indicating  that  a  rather 
low  order  quadrature  will  suffice.  Table  3  contains  these  results.  The  third  trail  was 
designed  to  test  the  multislab  aspect  of  the  F/v  algorithm.  For  this  case,  the  slab  of  A  =  1 
and  c  =  0.9  was  partitioned  into  6  slabs  each  with  c  =  0.9  and  of  widths: 


Slab  Number 

A 

1 

0.05 

2 

0.05 

3 

0.1 

4 

0.3 

5 

0.25 

6 

0.25 

The  choice  of  these  widths  was  dictated  by  the  results  of  [1]  which  contained  evaluations 
at  i  =  0.05,0.1,0.5,0.75,  and  1.0.  Three  inner  iterations  for  each  N  were  imposed  with 
Lm  =  20  and  e  =  10~5.  The  quadrature  order,  Lm3,  for  the  integrals  involving  boundary 
fluxes  was  varied  (starting  at  5  and  ending  at  25  with  intervals  of  5)  with  the  reflected 
flux,  at  x  ~  0.2,  and  the  transmitted  flux,  at  x  =  0.5,  displayed  m  Table  4.  In  comparison 
with  the  results  of  [1]  for  ,r  =  0.5,  the  algorithm  results  are  off  at  most  one  unit  in  the  fifth 
place  for  Lm3  >  15. 

The  final  trial  investigated  how  the  computational  time  and  accuracy  was  influenced  by  the 
number  of  inner  iterations  (LIT).  Table  5  shows  that  with  zero  inner  iteration  ( LIT = 1 ) , 
the  G  slab  case  did  not  converge  for  the  N=60  approximation.  With  1  or  more  iterations 
( LIT  > 2 ) .  the  fluxes  converge  for  N=54  with  no  improvement  as  LIT  increased  above  2. 
The  computational  time  on  a  CONVEX  mainframe  relative  to  zero  iterations  (C'PUo)  is  also 
given,  indicating  no  advantage  to  require  more  than  2  inner  iterations. 
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Table  1:  Reflected  and  Transmitted  Fluxes  for  Variation  of  e 


A  =  1.0,  c  =  0.9,  Lm  =  30,  IBR  =  4 


e  =  10-2 

and  N  =  14 

CO 

1 

o 

II 

and  N  =  14 

ft 

Reflected  Flux 

Transmitted  Flux 

Reflected  Flux 

Transmitted  Flux 

1.0 

2.10001E-01 

1.90264E-01 

2.10001E-01 

1.90264E-01 

0.9 

2.23885E-01 

2.00704E-01 

2.23885E-01 

2.00704E-01 

0.8 

2.39472E-01 

2.11883E-01 

2.39472E-01 

2.11883E-01 

0.7 

2.56973E-01 

2.23633E-01 

2.56973E-01 

2.23633E-01 

0.6 

2.76545E-01 

2.35537E-01 

2.76545E-01 

2.35537E-01 

0.5 

2.98151E-01 

2.46678E-01 

2.98151E-01 

2.46678E-01 

0.4 

3.21267E-01 

2.55202E-01 

3.21267E-01 

2.55202E  01 

0.3 

3.44242E-01 

2.57585E-01 

3.44242E-01 

2.57585E-01 

0.2 

3.63302E-01 

2.48316E-01 

3.63302E-01 

2.48316E-01 

0.1 

3.72679E-01 

2.24054E-01 

3.72679E-01 

2.24054E-01 

lc-10 

3.59370E-01 

1.86136E-01 

3.59370E-01 

1.86136E-01 

1 

o 

fH 

II 

and  N  =  16 

e  =  10-5 

and  N  =  22 

ft 

Reflected  Flux 

Transmitted  Flux 

Reflected  Flux 

Transmitted  Flux 

1.0 

2.10001E-01 

1.90264E-01 

2.10001E-01 

1.90265E-01 

0.9 

2.238S5E-01 

2.00704E-01 

2.23885E-01 

2.00704E-01 

o.s 

2.39472E-01 

2.11883E-01 

2.39472E-01 

2.11S83E-01 

0.7 

2.56974E-01 

2.23634E-01 

2.56974E-01 

2.23633E-01 

0.G 

2.76545E-01 

2.35536E-01 

2.76545E-01 

2.3553CE-01 

0.5 

2.98151E-01 

2.4G679E-01 

2.98151E-01 

2.46679E-01 

0.4 

3.212C7E-01 

2.55202E-01 

3.21267E-01 

2.55202E-01 

0.3 

3.44241E-01 

2.575S2E-01 

3.44241E-01 

2.57582E-01 

0.2 

3.63306E-01 

2.48327E-01 

3.83306E-01 

2.48328E-U1 

0.1 

3.72670E-01 

2.24032E-01 

3.72669E-01 

2.24033E-01 

IE-10 

3.59371E-01 

1.86138E-01 

3.59371E-01 

1.86138E-01 
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Table  2:  Reflected  and  Transmitted  Fluxes  for  e  Variation  for  Interior  Point 


x  =  0.5,  A  -  1.0,  c  =  0.9,  Lm  =  30,  IBR  =  4 


1.0 


0.9 


0.8 


0.7 


0.6 


0.5 


0.4 


0.3 


0.2 


0.1 


IE-10 


Reflected  Flux 


1.08052E-01 


1.17325E-01 


1.28289E-01 


1.41428E-01 


1.57408E-01 


1.77150E-01 


2.01882E-01 


2.32986E-01 


2.70683E-01 


3.08155E-01 


3.32748E-01 


Transmitted  Flux 


1.43159E-01 


1.55000E-01 


1.68882E-01 


1.85331E-01 


2.05035E-01 


2.28845E-01 


2.57648E-01 


2.91656E-01 


3.27322E-01 


3.47036E-01 


3.32748E-01 


Reflected  Flux 


1.08052E-01 


1.17325E-01 


1.28289E-01 


1.41428E-01 


1.57408E-01 


1.77150E-01 


2.01882E-01 


2.32986E-01 


2.70683E-01 


3.08155E-01 


3.32748E-01 


Transmitted  Flux 


1.43159E-01 


1.55000E-01 


1.68882E-01 


1.85331E-01 


2.05035E-01 


2.28845E-01 


2.57648E-01 


2.91656E-01 


3.27322E-01 


3.47036E-01 


3.32748E-01 


1.0 


0.9 


0.8 


0.7 


0.6 


0.5 


0.4 


0.3 


0.2 


0.1 


IE-10 


e  = 

1 

o 

T— » 

1  1 

Reflected  Flux 

Transmitted  Flux 

1.08052E-01 

1.43159E-01 

1.17325E-01 

1.55000E-01 

1.28289E-01 

1.6S882E-01 

1.41428E-01 

1.85331E-01 

1.57408E-01 

2.05035E-01 

1.77150E-01 

2.28845E-01 

2.01882E-01 

2.57648E-01 

2.32986E-01 

2.91657E-01 

2.70683E-01 

3.27320E-01 

3.08155E-01 

3.47040E-01 

3.32750E-01 

3.32750E-01 
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Table  3:  Reflected  and  Transmitted  Fluxes  from  Lm  Variation 


V 

Reflected  Flux 

Transmitted  Flux 

Reflected  Flux 

Transmitted  Flux 

1.0 

2.11696E-01 

1.92163E-01 

2.09973E-01 

1.90221E-01 

0.9 

2.25496E-01 

2.02572E-01 

2.23856E-01 

2.00656E-01 

0.8 

2.41032E-01 

2.13761E-01 

2.39441E-01 

2.11830E-01 

0.7 

2.58512E-01 

2.25564E-01 

2.56940E-01 

2.23573E-01 

0.6 

2.78096E-01 

2.37571E-01 

2.76509E-01 

2.35468E-01 

0.5 

2.99759E-01 

2.48878E-01 

2.98111E-01 

2.46601E-01 

0.4 

3.22988E-01 

2.57634E-01 

3.21222E-01 

2.551 12E-01 

0.3 

3.46144E-01 

2.60292E-01 

3.44190E-01 

2.57478E-01 

0.2 

3.65427E-01 

2.51228E-01 

3.63252E-01 

2.48219E-01 

0.1 

3.74S72E-01 

2.26683E-01 

3.72627E-01 

2.23955E-01 

IE-10 

3.61469E-01 

1.88096E-01 

3.59387E-01 

1.86136E-01 

Lm  =  8 

Ln 

i  =  16 

/‘ 

Reflected  Flux 

Transmitted  Flux 

Reflected  Flux 

Transmitted  Flux 

1.0 

2.10001E-01 

1.90264E-01 

2.10001E-01 

1.90265E-01 

0.9 

2.23885E-01 

2.00704E-01 

2.23885E-01 

2.00704E-01 

O.S 

2.39472E-01 

2.118S3E-01 

2.39472E-01 

2.118S3E-01 

0.7 

2.5C973E-01 

2.23G33E-01 

2.56974E-01 

2.23633E-01 

0.6 

2.7G545E-01 

2.35536E-01 

2.76545E-01 

2.35536E-01 

0.5 

2.98151E-01 

2.46678E-01 

2.98151E-01 

2.46679E-01 

0.4 

3.212G7E-01 

2.55201E-01 

3.21267E-01 

2.55202E-01 

0.3 

3.44240E-01 

2.57581E-01 

3.44241E-01 

2.57582E-01 

0.2 

3.63306E-01 

2.4S326E-01 

3.63306E-01 

2.48328E-01 

0.1 

3.72669E-01 

2.24032E-01 

2.24033E-01 

IE-10 

3.59374E-01 

1.8G139E-01 

3.59371E-01 

1.8G138E-01 

Table  4:  Reflected  and  Transmitted  Fluxes  from  Lm3  Variation 


Six  Slabs,  c  =  0.9,  Lm  =  20,  c  =  10-5,  IBR  =  4 


Lm3  =  5 

Lm3  =  10 

A* 

Reflected  Flux 
(. x  =  0.2) 

Transmitted  Flux 
( x  =  0.5) 

Reflected  Flux 
(x  =  0.2) 

Transmitted  Flux 
(x  =  0.5) 

1.0 

1.73095E-01 

1.43249E-01 

1.73230E-01 

1.43151E-01 

0.9 

1.85873E-01 

1.55093E-01 

1.86017E-01 

1.54991E-01 

0.8 

2.00532E-01 

1.68977E-01 

2.00687E-01 

1.68872E-01 

0.7 

2.17448E-01 

1.85428E-01 

2.17615E-01 

1.85322E-01 

0.G 

2.3704GE-01 

2.05131E-01 

2.37227E-01 

2.05025E-01 

0.5 

2.59741E-01 

2.28934E-01 

2.59935E-01 

2.28833E-01 

0.4 

2.85724E-01 

2.57721E-01 

2.85931E-01 

2.57635E-01 

0.3 

3.14366E-01 

2.91694E-01 

3.14577E-01 

2.91644E-01 

0.2 

3.42795E-01 

3.27284E-01 

3.4297SE-01 

3.27307E-01 

0.1 

3.G4943E-01 

3.46957E-01 

3.64926E-01 

3.47028E-01 

IE-10 

3.7778SE-01 

3.32759E-01 

3.77350E-01 

3.32744E-01 

Lm3  =  15 

Lm3  =  20 

/* 

Reflected  Flux 
(x-  =  0.2) 

Transmitted  Flux 
(x  =  0.5) 

Reflected  Flux 
(x  =  0.2) 

Transmitted  Flux 
(x  =  0.5) 

1.0 

1.73234E-01 

1.43159E-01 

1.73234E-01 

1.43159E-01 

0.9 

1.S6021E-01 

1.55000E-01 

1.86021E-01 

1.55000E-01 

0.8 

2.00691E-01 

1.68882E-01 

2.00691E-01 

1.68882E-01 

0.7 

2.17620E-01 

1.85331E-01 

2.17620E-01 

1.S5331E-01 

0.0 

2.37232E-01 

2.05035E-01 

2.37232E-01 

2.05035E-01 

0.5 

2.59941E-01 

2.28845E-01 

2.59941E-01 

2.28845E-01 

0.4 

2.S5937E-01. 

2.57G4SE-01 

2.85937E-01 

2.57647E-01 

0.3 

3.14583E-01 

2.91657E-01 

3.14583E-01 

2.91G57E-01 

0.2 

3.42984E-01 

3.27320E-01 

3.42984E-01 

3.27320E-01 

0.1 

3.G4932E-01 

3.47040E-01 

3.64932E-01 

3.47040E-01 

IE- 10 

3.7734GE-01 

3.3275  IE-01 

3.77347E-01 

3.32751E-01 
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Table  4:  Reflected  and  Transmitted  Fluxes  from  Lm3  Variation  (continued) 
Six  Slabs,  c  =  0.9,  Lm  =  20,  e  =  10~5,  IBR  =  4 


Lrn3 

=  25 

n 

Reflected  Flux 

Transmitted  Flux 

(x  =  0.2) 

(x  =  0.5) 

1.0 

1.73233E-01 

1.43159E-01 

0.9 

1.86021E-01 

1.55000 E-01 

0.8 

2.00691E-01 

1.68881E-01 

0.7 

2.17620E-01 

1.85331E-01 

0.6 

2.37232E-01 

2.05035E-01 

0.5 

2.59941  E-01 

2.28844E-01 

0.4 

2.85937E-01 

2.57647E-01 

0.3 

3.14583E-01 

2.91657E-01 

0.2 

3. 42984 E-01 

3.27320E-01 

0.1 

3.64932E-01 

3.47040E-01 

IE-10 

3.77346E-01 

3.32750E-01 
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Table  5:  Variation  of  CPU  Time  with  LIT 


LIT 

CPU  /  CPUo 

Comments 

1 

1.00 

Did  not  converge  for  N  =  60 

2 

1.83 

Converged  for  N  =  54 

3 

2.01 

Converged  for  N  =  54 

4 

2.24 

Converged  for  N  =  54 

V  Code  Operation 


This  section  provides  the  operational  details  for  the  MSLAB  code,  including 

•  a  detailed  subprogram  level  flow  chart  along  with  an  explanation  of  the  purpose  of 
each  subprogram, 

•  explanation  of  the  input  deck  and  output  using  three  examples  to  regenerate  data  of 
[1],  and 

•  detailed  programming  notes  to  ensure  proper  use. 

V.A  Flow  Chart 

A  list  of  the  subprogram  names  and  the  functions  performed  by  each  follows: 

Mui.TlSl.AH:  Main  driver  program  for  the  F/v  Method  in  a  heterogeneous  medium. 

Si:  Sets  boundary  conditions  and  interface  type  for  each  slab. 

IN'I'O:  Determines  integrals  over  boundary  fluxes  required  for  the  reflection  and  transmis¬ 
sion  coefficients. 

Post:  Calls  post-processing  routines  for  boundary  and  interior  fluxes. 

Pl’2:  Post- processing  routine  for  boundary  fluxes. 

Pl’2/.:  Post -processing  routine  for  interior  fluxes. 

Fn.\:  Determines  interior  fluxes. 

Soi.UTB:  Calls  for  the  matrix  inversion  of  the  F,v  collocation  equations  and  calculates  the 
reflection  and  transmission  coefficients. 

Mk.I.B:  Determines  matrix  elements  and  inhomogeneous  terms. 

INTH:  Evaluation  of  integrals  over  boundary  fluxes  for  inhomogeneous  terms  required  in 
collocation  equations  for  boundary  fluxes. 

INTR7.:  Evaluation  of  integrals  over  boundary  fluxes  for  inhomogeneous  terms  required  in 
collocation  equations  for  interior  fluxes. 

BASK:  Evaluation  of  basis  functions. 
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DBASE:  Used  when  p  =  v  (AM=V)  in  determination  of  matrix  element  Ba. 

BlSECl:  Performs  zero  search  via  Bi-Section  method. 

GauLeG:  Determination  of  Gauss- Legendre  abscissae  and  weights  from  [2]. 

ALAME:  Dispersion  relation. 

TEX!’:  Normal  FORTRAN  exponential  function  with  underflow  set  to  zero. 

Ma'I'INV:  Matrix  inversion  routine  that  uses  LU  decomposition. 

LuDcMP:  Performs  LU  decomposition  from  [2]. 

LuBkSU:  Uses  back  substitution  to  find  answer  from  [2], 

SFcN:  Function  S(x.y)  [eq.  (14h)]. 

CFc’N:  Function  C(x,y)  [eq.  (14g)]. 

MCoi’Y:  Copies  from  one  2-D  matrix  to  another. 

ICoPY:  Copies  from  integer  vector  to  another. 

All  subprograms  listed  above  are  Subroutines  except  SFcN.  CFcN,  and  Tl'.X I*  which  are 
Functions. 


11(1 


V 

BISEC1 

A  LAME 

1 

l  . 1 

I 


GAULEG 
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V.B  Input  Description 


The  input  deck  required  to  execute  this  program  is  explained  in  this  section.  Suggested 
values  and  limitations  are  noted  in  parenthesis.  The  following  text  is  also  contained  in  the 
source  code. 

LINE  1  NS  NUMBER  OF  SLABS  (<INS) 

LM3  GAUSS-LEGENDRE  INTEGRATION  ORDER  FOR  EVALUATION 
OF  INTEGRALS  OVER  BOUNDARY  FLUXES  (<ILM) 

LIT  NUMBER  OF  PASSES  THROUGH  SLABS  FOR  EACH  FN  APPROX.  (2) 

LINE  2  LI  INITIAL  N  OF  FN  APPROXIMATION 

L2  FINAL  N  OF  FN  APPROXIMATION  (<NP) 

L12D  INCREMENT  BETWEEN  APPROXIMATIONS  (2) 

LINE  3  MM1  NUMBER  OF  DIRECTION  EDIT  POINTS  (<IMM1) 

IG  -1  USE  G-L  ANGULAR  EDIT  POINTS  (ORDER  MM1 ) 

1  SPECIFY  UNIFORM  ANGULAR  EDIT  POINTS 

NOTE: IF  IG  IS  1  ENTER 
LINE  4  ANO  FIRST  EDIT  DIRECTION 

AN1  LAST  EDIT  DIRECTION 

LINE  5  ERR1  DESIRED  RELATIVE  ERROR  (1.0E-04) 

LINE  6  AMO  SOURCE  DIRECTION 

SO  INTENSITY  OF  BOUNDARY  SOURCE 

NOTE: ENTER  LINE  7  FOR  EACH  OF  THE  NS  SLABS 
LINE  7  WSR  SINGLE  SCATTER  ALBEDO 

DR  SLAB  WIDTH 

IBR  1  SHIFTED  LEGENDRE  POLYNOMIAL  BASIS  FCN3  (PL(2*X-1)) 

2  MONOMIAL  BASIS  FCNS  (X**(XAR*L)) 

3  SHIFTED  MONOMIAL  BASIS  FCNS  ( (2* (X**XAR) - 1 ) **L) 

4  SHIFTED  LEGENDRE  POLYNOMIAL  BASIS  FCNS 
(PL (U) ,U=2*(X**XAR)-1) 

LMR  G-L  QUADRATURE  ORDER  FOR  MATRIX  ELEMENTS  (<ILM) 

XAR  PARAMETER  FOR  BASIS  FCNS  OR  GAMMA  WHEN  IBR 
IS  3  OR  4  (0.75) 

SIR  EXPONENTIAL  VOLUME  SOURCE  STRENGTH 
ALP1R  EXPONENT  FOR  VOLUME  SOURCE 
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986  METHODS  DEUELOPHENT  FOR  ELECtfeOH  'TRffRSPffRt(U>  CAfsPAhP 
UB  RESEARCH  CENTER  BUFFALO  NY  B  D  GANAPOL  APR  92 
RL**-TR-92-48 


UNCLASSIFIED 


NL 


NXXR  NUMBER  OF  EDIT  POINTS  IN  SLAB  (<NPX) 


V.B.l  Sample  Problem  Input 

The  sample  problems  are  designed  to  reproduce  the  data  contained  in  Table  2  of  [1].  With 
the  MSLAB  program,  there  are  several  ways  to  generate  the  data.  Three  different  input 
decks  are  shown  and  explained  below. 

Input  deck  1: 


•  One  Slab  of  width  (A)  1.0. 

•  Evenly  spaced  interior  spatial  edit  points  at  a  spacing  of  0.05. 


I  16  3 

II  41  2 
10  1 
0.0  1.0 
1 .0E-05 
1.0  0.5 

0.9  1.0  4  30  0.75  0.0  0.0  20 


Input  deck  2: 


•  Three  Slabs  with  A’s  of  0.1,  0.4,  and  0.5. 

•  One,  Three,  and  One  interior  spatial  edit  points  for  the  slabs. 


3  16  3 
11  41  2 
10  1 
0.0  1.0 
1 .OE-05 
1.0  0.5 

0.9  0.1  4  30  0.75  0.0  0.0  2 
0.9  0.4  4  30  0.75  0.0  0.0  4 
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0.9  0.5  4  30  0.75  0.0  0.0  2 


Input  deck  3: 


•  Six  Slabs  with  A’s  of  0.05,  0.05,  0.1,  0.3,  0.25,  and  0.25. 

•  No  interior  spatial  edit  points  for  any  slab. 


6  16  3 
11  41  2 
10  1 


0.0  1.0 
1 .0E-05 
1.0  0.5 


0.9 

0.05 

4 

30 

0.75 

0.0 

0.0 

1 

0.9 

0.05 

4 

30 

0.75 

0.0 

0.0 

1 

0.9 

0.1 

4 

30 

0.75 

0.0 

0.0 

1 

0.9 

0.3 

4 

30 

0.75 

0.0 

0.0 

1 

0.9 

0.25 

4 

30 

0.75 

0.0 

0.0 

1 

0.9 

0.25 

4 

30 

0.75 

0.0 

0.0 

1 

V.C  Output  Description 


The  information  concerning  the  program’s  output  and  where  that  output  is  to  be  found 
is  contained  in  this  section.  It  also  contains  the  output  generated  by  the  example  input 
decks  of  the  preceding  section. 

The  output  files  are  as  found  in  Table  6.  The  contents  of  the  output  files  axe  found 
in  Table  7.  The  screen  output  (also  written  to  file  number  26  named  06)  shows  the 
convergence  criterion  during  the  program’s  execution.  The  first  column  of  screen  output 
(NN  =  )  contains  the  N  of  the  F/v  approximation.  The  second  column  (IE  =  )  shows  the 
number  of  boundary  points  that  have  converged.  The  third  column  (IEX  =  )  shows  the 
number  of  interior  points  that  have  converged.  The  final  entries  are  the  maximum  number 
of  boundary  and  interior  points  required  to  converge. 
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Table  6:  List  of  Files  Generated  by  MSLAB 


File  Number 

File  Name 

File  Contents 

21 

ol 

Transmission  and  Reflection  Coeff 

22 

o2 

Boundary  Fluxes  For  Each  Slab 

23 

o3 

Interior  Edit  (NXXR  points)  by  Slab 

26 

06 

Screen  Output 

Table  7:  Contents  of  Files  Generated  by  MSlab 


File 

Number 

File 

Name 

File  Contents  with  the  Text  Right  of  the  Arrow 

Describing  the  Data  in  that  Column 

21 

ol 

Slab  number  — ♦  Reflected  Coeff,  Transmitted  Coeff 

22 

o2 

Mu  vs  Boundary  Flux  — »  abs(Mu),  Flux  for  Mu<0,  Flux  for  Mu>0 

23 

o3 

Mu  vs  Interior  Flux  — ►  abs(Mu),  Flux  for  Mu<0,  Flux  for  Mu>0 

B34 


w*- 


V.C.l  Sample  Problem  Output 


The  program’s  output  due  to  the  input  decks  from  Section  V.B.l  are  shown  below  along 
with  plotted  results  that  match  Figure  1  in  Reference  [1]  (See  Figure  1  in  this  manual). 


Sample  Problem  1 


File  06: 


NN= 

11 

IE= 

0 

IEX= 

0 

22 

418 

NN= 

13 

IE= 

16 

IEX= 

0 

22 

418 

NN= 

15 

IE= 

18 

IEX= 

0 

22 

418 

NN= 

17 

IE= 

22 

IEX= 

0 

22 

418 

NN= 

19 

IE= 

22 

IEX= 

0 

22 

418 

NN= 

21 

IE= 

22 

IEX= 

-  404 

22 

418 

NN= 

23 

IE= 

22 

IEX= 

412 

22 

418 

NN= 

25 

IE= 

22 

IEX= 

416 

22 

418 

NN= 

27 

IE= 

22 

IEX= 

418 

22 

418 

File  ol: 


SLAB 

1 .337052E-01 
1 .337052E-01 
1 .337052E-01 
1 .337052E-01 
1 .337052E-01 
1 . 337052E-01 
1 . 337052E-01 
1 .337052E-01 
1 .337052E-01 

File  o2: 

FN  SOLUTION 
BDRY  FLUXES 


1 

2 . 958126E-01 
2 . 958125E-01 
2.958125E-01 
2.958125E-01 
2 . 958125E-01 
2.958125E-01 
2 . 958125E-01 
2 . 958125E-01 
2.958125E-01 
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SLAB  NR= 

1  X= [  O.OOOOOOOOOOOOOOOOE+OO 

1 . 00000E+00 

2 . 10001E-01 

1 . 90265E-01 

9.00000E-01 

2 . 23885E-01 

2 . 00704E-01 

8.00000E-01 

2.39472E-01 

2 . 11883E-01 

7.00000E-01 

2 . 56974E-01 

2 . 23633E-01 

6.00000E-01 

2 . 76545E-01 

2.35536E-01 

5.00000E-01 

2.98151E-01 

2 . 46679E-01 

4 . 00000E-01 

3.21267E-01 

2 . 55202E-01 

3.00000E-01 

3.44241E-01 

2 . 57582E-01 

2.00000E-01 

3 . 63306E-01 

2.48328E-01 

1 . 00000E-01 

3.72669E-01 

2 . 24033E-01 

1 .OOOOOE-IO 

3 . 5937 IE-01 

1 . 86138E-01 

File  o3: 

SLAB  NR= 

1  X=[  O.OOOOOOOOOOOOOOOOE+OO 

0 . OOOOOOOOOOOOOOOOE+OO 

- 

1 .OOOOOE+OO 

2. 10001E-01 

0. OOOOOE+OO 

9 . 00000E-01 

2 . 23885E-01 

0. OOOOOE+OO 

8 . 00000E-01 

2 . 39472E-01 

0. OOOOOE+OO 

7 . 00000E-01 

2 . 56974E-01 

0.  OOOOOE+OO 

6 . 00000E-01 

2.76545E-01 

0.  OOOOOE+OO 

5 . 00000E-01 

2.98151E-01 

0. OOOOOE+OO 

4 . 00000E-01 

3.21267E-01 

0. OOOOOE+OO 

3 . OOOOOE-Ol 

3.44241E-01 

0. OOOOOE+OO 

2 . 00000E-01 

3 . 63306E-01 

0. OOOOOE+OO 

1 .OOOOOE-Ol 

3.72669E-01 

0. OOOOOE+OO 

1 .OOOOOE-IO 

3 . 5937 IE-01 

0. OOOOOE+OO 

5 . OOOOOOOOOOOOOOOOE-02 

1 .OOOOOE+OO 

2.01846E-01 

1 . 80042E-02 

9. OOOOOE-Ol 

2 . 15593E-01 

1 .99500E-02 

8. OOOOOE-Ol 

2 . 311 16E-01 

2 . 23673E-02 

7.  OOOOOE-Ol 

2.48678E-01 

2 . 54507E-02 

6. OOOOOE-Ol 

2 . 68510E-01 

2 . 95197E-02 

5.  OOOOOE-Ol 

2.90702E-01 

3 . 51360E-02 

4. OOOOOE-Ol 

3 . 14917E-01 

4.33880E-02 

3. OOOOOE-Ol 

3 . 39767E-01 

5 . 66940E-02 

2.  OOOOOE-Ol 

3 . 61744E-01 

8 . 17099E-02 

1 .OOOOOE-Ol 

3 . 75369E-01 

1.45458E-01 

1  .OOOOOE-IO 

3 . 74858E-01 

3 . 74858E-01 

0.1000000000000000 


1.000000000000000  ] 


1.000000000000000  ] 
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l.OOOOOE+OO 

1 . 92839E-01 

3 . 55402E-02 

9.00000E-01 

2 . 06342E-01 

3 . 92759E-02 

8 . 00000E-01 

2 . 21674E-01 

4.38879E-02 

7.00000E-01 

2 . 39138E-01 

4.97251E-02 

6.00000E-01 

2 . 59037E-01 

5.73489E-02 

5 . 00000E-01 

2.81572E-01 

6 . 77242E-02 

4 . 00000E-01 

3.06584E-01 

8 . 26581E-02 

3 . 00000E-01 

3.32927E-01 

1 . 05960E-01 

2 . 00000E-01 

3 . 57282E-01 

1.47169E-01 

1 .00000E-01 

3 . 74075E-01 

2 . 36851E-01 

1 .00000E-10 

3.79338E-01 

3 . 79338E-01 

0.1500000000000000 

1.00000E+00 

1 . 83255E-01 

5 . 23288E-02 

9 . 00000E-01 

1.96434E-01 

5 . 76767E-02 

8.00000E-01 

2. 11477E-01 

6 . 42384E-02 

7 . 00000E-01 

2 . 28725E-01 

7.24779E-02 

6 . 00000E-01 

2.48544E-01 

8 . 31288E-02 

5 . 00000E-01 

2.71244E-01 

9 . 74200E-02 

4. OOOOOE-Ol 

2 . 96840E-01 

1 . 17571E-01 

3 . 00000E-01 

3 . 24431E-01 

1.47996E-01 

2.00000E-01 

3 . 50896E-01 

1 . 98623E-01 

1 .00000E-01 

3 . 70388E-01 

2 . 93091E-01 

1 .00000E-10 

3.79679E-01 

3 . 79679E-01 

0.2000000000000000 

1.00000E+00 

1 . 73234E-01 

6 . 82457E-02 

9 . 00000E-01 

1.86021E-01 

7 . 50244E-02 

8.00000E-01 

2 . 00691E-01 

8 . 32899E-02 

7 . 00000E-01 

2. 17620E-01 

9.35871E-02 

6 . 00000E-Q1 

2 . 37232E-01 

1 .06760E-01 

5 . 00000E-01 

2 . 59941E-01 

1 .24185E-01 

4.00000E-01 

2.85937E-01 

1 .48251E-01 

3 . 00000E-01 

3. 14583E-01 

1 .83409E-01 

2 . 00000E-01 

3 . 42984E-01 

2 . 38445E-01 

1 .OOOOOE-Ol 

3 . 64932E-01 

3 . 26737E-01 

1 .00000E-10 

3 . 77347E-01 

3 . 77347E-01 

0.2500000000000000 

1.00000E+00 

1 .62870E-01 

8 . 32219E-02 

9. 00000E-01 

1 . 75205E-01 

9 . 12525E-02 

8. 00000E-01 

1 .89426E-01 

1 .00983E-01 

7 .00000E-01 

2 . 05940E-01 

1. 13008E-01 

6. OOOOOE-Ol 

2 . 25227E-01 

1 .28232E-01 

5.  OOOOOE-Ol 

2.47801E-01 

1 . 48082E-01 
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4 . 00000E-01 

2 . 74027E-01 

1 . 74930E-01 

3 . 00000E-01 

3 . 03554E-01 

2. 12864E-01 

2 . 00000E-01 

3 . 33770E-01 

2 . 68706E-01 

1 . 00000E-01 

3 . 58066E-01 

3.45790E-01 

1 .00000E-10 

3 . 73044E-01 

3.73044E-01 

0.3000000000000000 

1 .00000E+00 

1 . 52238E-01 

9.72184E-02 

9 . 00000E-01 

1 . 64062E-01 

1 . 06327E-01 

8 . 00000E-01 

1 . 77763E-01 

1 . 17294E-01 

7 . 00000E-01 

1 . 93773E-01 

1 . 30738E-01 

6 . OOOOOE-Ol 

2 . 12622E-01 

1 .47578E-01 

5.00000E-01 

2.34919E-01 

1 . 69218E-01 

4. OOOOOE-Ol 

2 . 61209E-01 

1.97871E-01 

3. 00000E-01 

2 . 91446E-01 

2 . 37011E-01 

2. OOOOOE-Ol 

3 . 23378E-01 

2.91137E-01 

1 .OOOOOE-Ol 

3 . 50015E-01 

3 . 55312E-01 

1 .00000E-10 

3 . 67190E-01 

3 . 67190E-01 

0.3500000000000000 

1 .OOOOOE+OO 

1 . 41393E-01 

1 . 10215E-01 

9. OOOOOE-Ol 

1 .52653E-01 

1 .20235E-01 

8. OOOOOE-Ol 

1 . 65767E-01 

1 . 32222E-01 

7. OOOOOE-Ol 

1 . 81 186E-01 

1 .46797E-01 

6. OOOOOE-Ol 

1.99486E-01 

1 . 64856E-01 

5. OOOOOE-Ol 

2 . 21366E-01 

1.87722E-01 

4. OOOOOE-Ol 

2 . 47550E-01 

2.17351E-01 

3. OOOOOE-Ol 

2 . 78320E-01 

2 . 56448E-01 

2. OOOOOE-Ol 

3 . 11877E-01 

3 . 07160E-01 

1 .OOOOOE-Ol 

3 . 40930E-01 

3.58505E-01 

1 .00000E-10 

3 . 60063E-01 

3 . 60063E-01 

0.4000000000000000 

1 .OOOOOE+OO 

1 . 30386E-01 

1 . 22202E-01 

9 .OOOOOE-Ol 

1 . 41032E-01 

1 .32976E-01 

8. OOOOOE-Ol 

1 . 53493E-0 1 

1.45780E-01 

7. OOOOOE-Ol 

1 . 68236E-01 

1 .61218E-01 

6.  OOOOOE-Ol 

1.85876E-01 

1.80138E-01 

5. OOOOOE-Ol 

2.07195E-01 

2.03734E-01 

4. OOOOOE-Ol 

2 . 33093E-01 

2.33637E-01 

3.  OOOOOE-Ol 

2 . 64203E-01 

2.71721E-01 

2.  OOOOOE-Ol 

2 . 99284E-01 

3.17935E-01 

1  .OOOOOE-Ol 

3 . 30905E-01 

3 . 57403E-01 

1 .00000E-10 

3 . 51858E-01 

3.51858E-01 

0.4500000000000000 
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l.OOOOOE+OO 

1 . 19259E-01 

1.33181E-01 

9 . 00000E-01 

1 . 29242E-01 

1.44560E-01 

8 . 00000E-01 

1.40986E-01 

1 . 57991E-01 

7 . 00000E-01 

1 . 54969E-01 

1.74047E-01 

6 . 00000E-01 

1 . 7 1837E-01 

1 . 93504E-01 

5 . 00000E-01 

1 .92446E-01 

2 . 17394E-01 

4 . 00000E-01 

2 . 17865E-01 

2.46989E-01 

3.00000E-01 

2.49098E-01 

2 . 83313E-01 

2 . 00000E-01 

2.85575E-01 

3 . 24402E-01 

1 .00000E-01 

3 . 19983E-01 

3 . 53302E-01 

1 .00000E-10 

3 . 42719E-01 

3.42719E-01 

1.5000000000000000 

1.00000E+00 

1 . 08052E-01 

1.43159E-01 

9 . 00000E-01 

1 . 17325E-01 

1 .55000E-01 

8.00000E-01 

1 .28289E-01 

1 .68882E-01 

7 . 00000E-01 

1 .41428E-01 

1 . 85331E-01 

6 . OOOOOE-Ol 

1 .57408E-01 

2 . 05035E-01 

5 .OOOOOE-Ol 

1.77150E-01 

2 . 28845E-01 

4 . OOOOOE-Ol 

2 . 01882E-01 

2 . 57648E-01 

3 . 00000E-01 

2 . 32986E-01 

2 . 91657E-01 

2. OOOOOE-Ol 

2 . 70683E-01 

3 . 27320E-01 

1. OOOOOE-Ol 

3.08157E-01 

3.47040E-01 

1  .OOOOOE-IO 

3 . 32750E-01 

3 . 32750E-01 

1.5500000000000000 

1 . OOOOOE+OO 

9 . 6801 IE-02 

1 . 52145E-01 

9. OOOOOE-Ol 

1.05319E-01 

1.64317E-01 

8. OOOOOE-Ol 

1.  15441E-01 

1 .78485E-01 

7. OOOOOE-Ol 

1.27650E-01 

1 .95124E-01 

6. OOOOOE-Ol 

1 . 42624E-01 

2. 14817E-01 

5. OOOOOE-Ol 

1 .61333E-01 

2.38220E-01 

4.  OOOOOE-Ol 

1 .85148E-01 

2 . 65836E-01 

3.  OOOOOE-Ol 

2 . 15827E-01 

2 . 97129E-01 

2. OOOOOE-Ol 

2 . 54496E-01 

3 . 27299E-01 

1 .OOOOOE-Ol 

2 . 95353E-01 

3 . 39157E-01 

1 .OOOOOE-IO 

3 . 22030E-01 

3 . 22030E-01 

1.6000000000000000 

l.OOOOOE+OO 

8.55408E-02 

1 . 60152E-01 

9. OOOOOE-Ol 

9 . 32583E-02 

1 . 72531E-01 

8. OOOOOE-Ol 

1 .02478E-01 

1 .86836E-01 

7. OOOOOE-Ol 

1  .  13671E-01 

2 . 03478E-01 

6 .OOOOOE-Ol 

1 .27517E-01 

2 . 22931E-01 

5 .OOOOOE-Ol 

1 .45017E-01 

2.45648E-01 
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4 . 00000E-01 

1.67660E-01 

2.71761E-01 

3 . 00000E-01 

1 . 97563E-01 

3 . 00059E-01 

2 . 00000E-01 

2 . 36853E-01 

3 . 24831E-01 

1 .00000E-01 

2 . 81407E-01 

3 . 30007E-01 

1 . 00000E-10 

3. 10610E-01 

3. 10610E-01 

'.6500000000000000 

1.00000E+00 

7.43058E-02 

1 . 67195E-01 

9 . 00000E-01 

8.11804E-02 

1.79666E-01 

8 . 00000E-01 

8 . 94368E-02 

1.93969E-01 

7 . 00000E-01 

9 . 95276E-02 

2 . 10446E-01 

6 . 00000E-01 

1 . 12119E-01 

2.29456E-01 

5 . 00000E-01 

1 .28221E-01 

2 . 51251E-01 

4 . 00000E-01 

1 .49407E-01 

2 . 75607E-01 

3 . 00000E-01 

1 . 781 16E-01 

3 . 00734E-01 

2 . 00000E-01 

2. 17534E-01 

3 . 20305E-01 

1 .OOOOOE-Ol 

2 . 66023E-01 

3. 19821E-01 

1 .00000E-10 

2 . 98522E-01 

2 . 98522E-01 

1.7000000000000000 

1 .00000E+00 

6.31312E-02 

1.73289E-01 

9. 00000E-01 

6 . 91218E-02 

1 . 85745E-01 

8.  OOOOOE-Ol 

7 . 63552E-02 

1.J9916E-01 

7. OOOOOE-Ol 

8 . 52560E-02 

2 . 16077E-01 

6. OOOOOE-Ol 

9 . 64610E-02 

2.34466E-01 

5. OOOOOE-Ol 

1 . 10963E-01 

2.55138E-01 

4. OOOOOE-Ol 

1 . 30376E-01 

2 . 77542E-01 

3.  OOOOOE-Ol 

1 . 57393E-01 

2 . 99397E-01 

2. OOOOOE-Ol 

1 .96250E-01 

3 . 14030E-01 

1  .OOOOOE-Ol 

2 . 48696E-01 

3 . 08746E-01 

1 .00000E-10 

2.85774E-01 

2.85774E-01 

:  .7500000000000000 

1  .OOOOOE+OO 

5 . 20545E-02 

1 .78448E-01 

9. OOOOOE-Ol 

5.71217E-02 

1 .90787E-01 

8. OOOOOE-Ol 

6 . 32737E-02 

2 . 047 1 IE-01 

7. OOOOOE-Ol. 

7 . 08963E-02 

2 . 20417E-01 

6.  OOOOOE-Ol 

8 . 05794E-02 

2 . 38028E-01 

5.  OOOOOE-Ol 

9 . 32667E-02 

2 . 57409E-01 

4.  OOOOOE-Ol 

1 . 10552E-01 

2 . 777 1  IE-01 

3.  OOOOOE-Ol 

1 . 35283E-01 

2 . 96255E-01 

2. OOOOOE-Ol 

1 . 72635E-01 

3 . 06245E-01 

1. OOOOOE-Ol 

2 . 28599E-01 

2 . 96869E-01 

1 .00000E-10 

2 . 72344E-01 

2 . 72344E-01 

0.8000000000000000 
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1 .OOOOOE+OO 

4.11169E-02 

1.82682E-01 

9.00000E-01 

4 . 52238E-02 

1.94812E-01 

8 . 00000E-01 

5 . 02378E-02 

2 . 08379E-01 

7 . 00000E-01 

5 . 64944E-02 

2 . 23507E-01 

6 . 00000E-01 

6 . 45157E-02 

2 . 40203E-01 

5 .00000E-01 

7.51588E-02 

2 . 58151E-01 

4 . 00000E-01 

8 . 99252E-02 

2.76238E-01 

3 . 00000E-01 

1.11663E-01 

2 . 91477E-01 

2 . 00000E-01 

1.46225E-01 

2 . 97128E-01 

1 .00000E-01 

2 . 04393E-01 

2 . 84226E-01 

1 .00000E-10 

2 . 58174E-01 

2 . 58174E-01 

.8500000000000000 

1 .00000E+00 

3 . 03663E-02 

1 .85998E-01 

9 . 00000E-01 

3 . 34794E-02 

1 . 97830E-01 

8 . 00000E-01 

3 . 73017E-02 

2. 10940E-01 

7 . 00000E-01 

4 . 2 1058E-02 

2 . 25378E-01 

6 . 00000E-01 

4 . 83234E-02 

2 . 41037E-01 

5 . OOOOOE-Ol 

5 . 6680  IE-02 

2 . 57433E-01 

4 . 00000E-01 

6 . 84927E-02 

2 . 73223E-01 

3. OOOOOE-Ol 

8 . 6403  IE-02 

2 . 85190E-01 

2 .OOOOOE-Ol 

1 . 16452E-01 

2 . 86797E-01 

1 .OOOOOE-Ol 

1  73925E-01 

2 . 70803E-01 

1 .00000E-10 

2 . 43137E-01 

2.43137E-01 

1. 9000000000000000 

1 .OOOOOE+OO 

1 .98624E-02 

1 .88393E-01 

9. OOOOOE-Ol 

2. 19531E-02 

1 .99844E-01 

8. OOOOOE-Ol 

2 . 45350E-02 

2. 12403E-01 

7. OOOOOE-Ol 

2 . 78041E-02 

2 . 26046E-01 

6. OOOOOE-Ol 

3 . 20761E-02 

2 . 40559E-01 

5. OOOOOE-Ol 

3 . 78942E-02 

2 . 55302E-01 

4. OOOOOE-Ol 

4 . 62777E-02 

2 . 68731E-01 

3. OOOOOE-Ol 

5 . 93824E-02 

2 . 77476E-01 

2. OOOOOE-Ol 

8 . 26387E-02 

2 . 75302E-01 

1 .OOOOOE-Ol 

1 .33762E-01 

2 . 56518E-01 

1 .00000E-10 

2 . 26978E-01 

2 . 26978E-0 1 

.9500000000000000 

1 .OOOOOE+OO 

9.68862E-03 

1 .89844E-01 

9. OOOOOE-Ol 

1 .07361E-02 

2 . 00832E-01 

8. OOOOOE-Ol 

1 .20374E-02 

2. 12750E-01 

7. OOOOOE-Ol 

1 .36976E-02 

2 . 25499E-01 

6. OOOOOE-Ol 

1 . 58887E-02 

2 . 38763E-01 

5. OOOOOE-Ol 

1 .89135E-02 

2 . 5 1759E-01 
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4 . 00000E-01 

2 . 3359  IE-02 

2.62774E-01 

3 . 00000E-01 

3.05306E  -02 

2 . 68341E-01 

2 . 00000E-01 

4.40245E-02 

2 . 62594E-01 

1 .00000E-01 

7 . 84908E-02 

2 . 41 157E-0 1 

1 .00000E-10 

2 . 09079E-01 

2 . 09079E-01 

1.000000000000000 

1 .00000E+00 

0 . OOOOOE+OO 

1 .90265E-01 

9 . 00000E-01 

0 . OOOOOE+OO 

2.00704E-01 

8.00000E-01 

0. OOOOOE+OO 

2.  U883E-01 

7 . OOOOOE-Ol 

0. OOOOOE+OO 

2 . 23633E-01 

6 . OOOOOE-Ol 

0.  OOOOOE+OO 

2 . 35536E-01 

5. OOOOOE-Ol 

0.  OOOOOE+OO 

2.46679E-01 

4. OOOOOE-Ol 

0.  OOOOOE+OO 

2 . 55202E-01 

3. OOOOOE-Ol 

0.  OOOOOE+OO 

2 . 57582E-01 

2. OOOOOE-Ol 

0.  OOOOOE+OO 

2 . 48328E-01 

1 .OOOOOE-Ol 

0.  OOOOOE+OO 

2 . 24033E-01 

1  .  OOOOOE- 10 

0.  OOOOOE+OO 

1 .86138E-01 
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Sample  Problem  2 


File  c>6: 


N»= 

u 

IE= 

0  IEX= 

0 

66 

110 

NN* 

ii 

IE= 

0  IEX= 

0 

66 

110 

NN= 

ii 

IE= 

0  IEX= 

0 

66 

110 

NN= 

13 

IE= 

0  IEX= 

0 

66 

110 

NN  = 

13 

IE= 

0  IEX= 

0 

66 

110 

NN  = 

13 

IE= 

0  IEX= 

0 

66 

110 

NN= 

15 

IE= 

14  IEX= 

0 

66 

110 

NN= 

15 

IE= 

15  IEX= 

0 

66 

110 

NN  = 

15 

IE= 

15  IEX= 

0 

66 

110 

NN= 

17 

IE= 

51  IEX= 

0 

66 

110 

NN  = 

17 

IE= 

50  IEX= 

0 

66 

110 

NN  = 

17 

IE= 

50  IEX= 

0 

66 

110 

NN= 

19 

IE= 

59  IEX= 

0 

66 

110 

NN= 

19 

IE= 

59  IEX= 

0 

66 

no 

NN= 

19 

IE= 

5?  IEX= 

0 

66 

no 

NN  = 

21 

IE= 

63  IEX= 

0 

66 

no 

NN  = 

21 

IE= 

63  IEX= 

0 

66 

no 

NN= 

21 

IE= 

63  IEX= 

0 

66 

no 

NN  = 

23 

IE= 

66  IEX= 

0 

66 

no 

NN= 

23 

IE= 

66  IEX= 

0 

66 

no 

NN= 

25 

IE= 

66  IEX= 

109 

66 

no 

=2 

z: 

u 

27 

IE= 

66  IEX= 

110 

66 

no 

File  <>1: 


SLAB 

1 . 285502E-01 
1 .336903E-01 
1  337049E-01 
1 .337051E-01 
1 .337052E-01 
1 .337052E-01 
1 .337052E-01 
1 .337052E-01 
1 .337052E-01 


1 

4 . 832700E-0 1 
4 . 837045E-01 
4 . 837057E-01 
4 . 837056E-01 
4 . 837056E-01 
4 . 837056E-01 
4.837056E-01 
4 . 837056E-01 
4 . 837056E-01 
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V.D  Detailed  Program  Notes 

rV<at  oilorl  nmirrarmnimr  nnf pc  are  on  vnn  in  liotps  rpfpr  In  rntiminiifc  nti/1 


SLAB 

1 . 248301E-01 
1 . 257041E-01 
1 . 257066E-01 
1 . 257069E-01 
1 . 257070E-01 
1 . 257070E-01 
1 . 257070E-01 
1 . 257069E-01 
1 . 257069E-01 

SLAB 

7 . 809918E-02 
7 . 824439E-02 
7 . 824487E-02 
7 . 824500E-02 
7 . 824504E-02 
7.824503E-02 
7 . 824503E-02 
7 . 824502E-02 
7 . 824502E-02 

File  o2: 


2 

4 . 033682E-01 
4 . 039195E-01 
4 . 039214E-01 
4 . 039218E-01 
4 . 039219E-01 
4 . 039219E-01 
4 . 039219E-01 
4 . 039218E-01 
4 . 039218E-01 

3 

2 . 954598E-01 
2 . 958111E-01 
2 . 958122E-01 
2 . 958125E-01 
2.958126E-01 
2 . 958126E-01 
2 . 958126E-01 
2.958126E-01 
2 . 958126E-01 


FN  SOLUTION 

BDRY  FLUXES 

SLAB  NR= 

i  x=[  o.oooooooooooooooo; 

1 .OOOOOE+OO 

2. 10001E-01 

3 . 55403E-02 

9 . 00000E-01 

2 . 23885E-01 

3 . 92760E-02 

8 . OOOOOE-Ol 

2 . 39472E-01 

4 . 38880E-02 

7 . 00000E-01 

2 . 56974E-01 

4 . 97252E-02 

6. OOOOOE-Ol 

2 . 76545E-01 

5 . 73490E-02 

5 .OOOOOE-Ol 

2 . 98152E-01 

6 . 77244E-02 

4. OOOOOE-Ol 

3.21267E-01 

8 . 26583E-02 

3. OOOOOE-Ol 

3.44241E-01 

1 .05961E-01 

2. OOOOOE-Ol 

3 . 63307E-01 

1 .47170E-01 

1 .OOOOOE-Ol 

3.72670E-01 

2.36851E-01 

1 .OOOOOE-IO 

3.59367E-01 

3 . 79340E-01 

SLAB  NR= 

2  X=[  0.1000000000000000 

1 .OOOOOE+OO 

1 .92839E-01 

1 .43159E-01 

0.1000000000000000 


0.5000000000000000 


,  l.J  ; 

9.00000E-01 

2.06342E-01 

1.55000E-01 

8.00000E-01 

2 . 21674E-01 

1.68882E-01 

7.00000E-01 

2 . 39138E-01 

1.85332E-01 

6 . 00000E-01 

2 . 59037E-01 

2 . 05036E-01 

5 . 00000E-01 

2 . 81572E-01 

2 . 28845E-01 

4.00000E-01 

3 . 06585E-01 

2 . 57648E-01 

3 .  OOOOOE-Ol 

3 . 32927E-01 

2 . 91657E-01 

2 . 00000E-01 

3 . 57282E-01 

3.27320E-01 

1.  OOOOOE-Ol 

3.74076E-01 

3 . 47040E-01 

1  .OOOOOE-IO 

3 . 79344E-01 

3.32751E-01 

SLAB  NR= 

3  X= [  0.5000000000000000 

1.000000000000000 

] 

1  .OOOOOE+OO 

1.08052E-01 

1 . 90265E-01 

9.  OOOOOE-Ol 

1 . 17325E-01 

2.00704E-01 

8.  OOOOOE-Ol 

1 . 28289E-01 

2.11883E-01 

7.  OOOOOE-Ol 

1.41428E-01 

2 . 23634E-01 

6.  OOOOOE-Ol 

1 . 57408E-01 

2 . 35536E-01 

5.  OOOOOE-Ol 

1.77150E-01 

2.46679E-01 

4.  OOOOOE-Ol 

2 . 01882E-01 

2.55202E-01 

3.  OOOOOE-Ol 

2 . 32986E-01 

2.57582E-01 

2.  OOOOOE-Ol 

2 . 70683E-01 

2.48328E-01 

1 .OOOOOE-Ol 

3 . 08157E-01 

2.24033E-01 

1 .OOOOOE-IO 

3 . 32751E-01 

1 .86138E-01 

File  o3: 

SLAB  NR= 

1  X=[  0. OOOOOOOOOOOOOOOOE+OO, 

0.1000000000000000 

] 

0 . OOOOOOOOOOOOOOOOE+OO 

1.  OOOOOE+OO 

2. 10001E-01 

0.  OOOOOE+OO 

9.  OOOOOE-Ol 

2 . 23885E-01 

0.  OOOOOE+OO 

8.  OOOOOE-Ol 

2 . 39472E-01 

0.  OOOOOE+OO 

7. OOOOOE-Ol 

2 . 56974E-01 

0.  OOOOOE+OO 

6.  OOOOOE-Ol 

2 . 76545E-01 

0.  OOOOOE+OO 

5. OOOOOE-Ol 

2 . 98152E-01 

0. OOOOOE+OO 

4.  OOOOOE-Ol 

3 . 21267E-01 

0.  OOOOOE+OO 

3. OOOOOE-Ol 

3.44241E-01 

0.  OOOOOE+OO 

2.  OOOOOE-Ol 

3 . 63307E-01 

0.  OOOOOE+OO 

1.  OOOOOE-Ol 

3 . 72670E-01 

0. OOOOOE+OO 

1.  OOOOOE-IO 

3 . 59367E-01 

0. OOOOOE+OO 

5 . 0000000000000000E-02 

1 .OOOOOE+OO 

2 . 01847E-01 

1 .80042E-02 

9. OOOOOE-Ol 

2 . 15593E-01 

1 .99500E-02 
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8 . 00000E-01 

2.31116E-01 

2.23673E-02 

7 . 00000E-01 

2.48678E-01 

2 . 54507E-02 

6 . 00000E-01 

2 . 68510E-01 

2 . 95197E-02 

5 . 00000E-01 

2 . 90702E-01 

3 . 51360E-02 

4.00000E-01 

3 . 14918E-01 

4.33880E-02 

3.00000E-01 

3 . 39768E-01 

5.66940E-02 

2 . 00000E-01 

3 . 61745E-01 

8.17100E-02 

1 .00000E-01 

3 . 75369E-01 

1.45458E-01 

1 -00000E-10 

3 . 74876E-01 

3.74876E-01 

SLAB  NR= 

2  X=[  0.1000000000000000 

,  0.5000000000000000 

0.1000000000000000 


1 .00000E+00 

1 . 92839E-01 

3.55403E-02 

9 . OOOOOE-Ol 

2 . 06342E-01 

3 . 92760E-02 

8.00000E-01 

2 . 21674E-01 

4.38880E-02 

7. 00000E-01 

2 . 39138E-01 

4 . 97252E-02 

6. OOOOOE-Ol 

2.59037E-01 

5.73490E-02 

5. OOOOOE-Ol 

2 . 81572E-01 

6.77243E-02 

4. OOOOOE-Ol 

3 . 06585E-01 

8 . 26583E-02 

3. OOOOOE-Ol 

3 . 32927E-01 

1 . 05960E-01 

2. OOOOOE-Ol 

3 . 57282E-01 

1.47170E-01 

1. OOOOOE-Ol 

3 . 74076E-01 

2.36851E-01 

1 .00000E-10 

3.79344E-01 

3.79339E-01 

0.2000000000000000 

1 . OOOOOE+OO  1 . 73234E-01 

6.82458E-02 

9. OOOOOE-Ol 

1 .86021E-01 

7.50246E-02 

8. OOOOOE-Ol 

2.00691E-01 

8.32901E-02 

7. OOOOOE-Ol 

2. 17620E-01 

9.35873E-02 

6. OOOOOE-Ol 

2 . 37232E-01 

1 .06760E-01 

5. OOOOOE-Ol 

2 . 59941E-01 

1 . 24186E-01 

4. OOOOOE-Ol 

2 . 85937E-01 

1.48252E-01 

3. OOOOOE-Ol 

3. 14583E-01 

1.83409E-01 

2. OOOOOE-Ol 

3.42984E-01 

2 . 38446E-01 

1 .OOOOOE-Ol 

3 . 64932E-01 

3.26739E-01 

1 .00000E-10 

3.77340E-01 

3.77340E-01 

0.3000000000000000 

1 . OOOOOE+OO  1 . 52238E-01 

9 . 72185E-02 

9. OOOOOE-Ol 

1 . 64062E-01 

1 . 06327E-01 

8. OOOOOE-Ol 

1 . 77763E-01 

1 . 17294E-01 

7. OOOOOE-Ol 

1 . 93773E-01 

1.30739E-01 

6.  OOOOOE-Ol 

2. 12622E-01 

1.47578E-01 

5.  OOOOOE-Ol 

2 . 34919E-01 

1 . 69218E-01 

4. OOOOOE-Ol 

2 . 61209E-01 

1.97871E-01 
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3.00000E-01 

2 . 91446E-01 

2.37011E-01 

2.00000E-01 

3 . 23379E-01 

2 . 91137E-01 

1.00000E-01 

3 . 50015E-01 

3 . 55312E-01 

1.00000E-10 

3.67192E-01 

3.67192E-01 

1.4000000000000000 

1.00000E+00 

1 . 30386E-01 

1 . 22202E-01 

9 . 00000E-01 

1.41032E-01 

1.32977E-01 

8 . 00000E-01 

1 . 53493E-01 

1.45780E-01 

7.00000E-01 

1.68236E-01 

1 . 61218E-01 

6.00000E-01 

1 . 85876E-01 

1.80138E-01 

5 . 00000E-01 

2 . 07195E-01 

2.03734E-01 

4 . 00000E-01 

2 . 33093E-01 

2 . 33638E-01 

3 . OOOOOE-Ol 

2 . 64203E-01 

2.71721E-01 

2 . 00000E-01 

2 . 99284E-01 

3. 17935E-01 

1 .00000E-01 

3 . 30906E-01 

3.57402E-01 

1 .00000E-10 

3 . 51854E-01 

3 . 51854E-01 

,AB  NR= 

3  X= [  0.5000000000000000 

.5000000000000000 

1 .OOOOOE+OO 

1.08052E-01 

1.43159E-01 

9. OOOOOE-Ol 

1 . 17325E-01 

1 . 55000E-01 

8. OOOOOE-Ol 

1 . 28289E-01 

1 .68882E-01 

7. OOOOOE-Ol 

1 . 41428E-01 

1 .85332E-01 

6. OOOOOE-Ol 

1 . 57408E-01 

2 . 05035E-01 

5. OOOOOE-Ol 

1 . 77150E-01 

2 . 28845E-01 

4. OOOOOE-Ol 

2.01882E-01 

2 . 57648E-01 

3. OOOOOE-Ol 

2.32986E-01 

2 . 91657E-01 

2. OOOOOE-Ol 

2.70683E-01 

3 . 27320E-01 

1 .OOOOOE-Ol 

3.08157E-01 

3.47040E-01 

1 .00000E-10 

3 . 32751E-01 

3 . 32751E-01 

7500000000000000 

1 .OOOOOE+OO 

5 . 20545E-02 

1.78448E-01 

9. OOOOOE-Ol 

5.71217E-02 

1 .90787E-01 

8. OOOOOE-Ol 

6.32737E-02 

2.04711E-01 

7. OOOOOE-Ol 

7 . 08963E-02 

2.20417E-01 

6. OOOOOE-Ol 

8 . 05794E-02 

2 . 38029E-01 

5. OOOOOE-Ol 

9 . 32667E-02 

2 . 57409E-01 

4. OOOOOE-Ol 

1 . 10552E-01 

2 . 7771 IE-01 

3. OOOOOE-Ol 

1 . 35283E-01 

2 . 96255E-01 

2. OOOOOE-Ol 

1 . 72635E-01 

3.06246E-01 

1. OOOOOE-Ol 

2 . 28599E-01 

2 . 96869E-01 

1 .00000E-10 

2 . 72345E-01 

2.72345E-01 

1.000000000000000 


1.000000000000000  ] 
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1 .OOOOOE+OO 
9 . 00000E-01 
8 . 00000E-01 
7 . 00000E-01 
6 . 00000E-01 
5 . 00000E-01 
4 . 00000E-01 
3 . 00000E-01 
2 . 00000E-01 
1 .OOOOOE-Ol 
1 .OOOOOE-IO 


0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 


1 . 90265E-01 
2 . 00704E-01 
2 . 11883E-01 
2 . 23634E-01 
2 . 35536E-01 
2.46679E-01 
2 . 55202E-01 
2 . 57582E-01 
2.48328E-01 
2 . 24033E-01 
1 . 86138E-01 


Sample  Problem  3 


File  06: 


NN= 

11  IE® 

0  I  EX® 

NN= 

11  IE® 

0  I  EX® 

NN® 

11  IE® 

0  I  EX® 

NN= 

13  IE® 

0  I  EX® 

NN= 

13  IE® 

0  I EX® 

NN* 

13  IE® 

0  I  EX® 

NN® 

15  IE® 

0  I EX® 

NN® 

15  IE® 

0  I EX® 

NN® 

15  IE® 

0  I  EX® 

NN® 

17  IE® 

0  I  EX® 

NN= 

17  IE® 

0  IEX® 

NN® 

17  IE® 

0  IEX® 

NN® 

19  IE® 

58  IEX® 

NN® 

19  IE® 

47  IEX® 

NN= 

19  IE® 

40  IEX® 

NN® 

21  IE® 

117  IEX* 

NN® 

21  IE® 

108  IEX* 

NN® 

21  IE® 

103  IEX* 

NN® 

23  IE® 

122  IEX* 

NN® 

23  IE® 

121  IEX* 

NN= 

23  IE® 

122  IEX* 

NN® 

25  IE® 

130  IEX* 

NN® 

25  IE® 

129  IEX* 

NN® 

25  IE® 

129  IEX* 

NN® 

27  IE® 

132  IEX* 

NN® 

27  IE® 

132  IEX* 

File  ol: 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 

132 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


SLAB 

3 . 937156E-02 
1 . 281899E-01 
1 . 333497E-01 
1 . 336872E-01 
1 . 337041E-01 


1 

4 . 881898E-01 
4 . 920301E-01 
4 . 922490E-01 
4 . 922631E-01 
4 . 922638E-01 
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1 . 337052E-01 
1 . 337054E-01 
1 . 337055E-01 
1 . 337056E-01 

SLAB 

7 . 580211E-02 
1 . 276844E-01 
1 . 299296E-01 
1 . 300631E-01 
1 . 300698E-01 
1 . 300703E-01 
1 . 300704E-01 
1 . 300704E-01 
1 . 300704E-01 

SLAB 

9 . 772624E-02 
1.246937E-01 
1 .256555E-01 
1 . 257041E-01 
1 . 257067E-01 
1.257069E-01 
1 . 257070E-01 
1  257071E-01 
1 . 25707 IE-01 

SLAB 

1 . 092035E-01 
1 . 151586E-01 
1 . 155566E-01 
1. 155764E-01 
1. 155775E-01 
1. 155776E-01 
1. 155777E-01 
1. 155777E-01 
1 . 155777E-01 

SLAB 

7 . 555241E-02 
7.807422E-02 
7 . 823655E-02 


4 . 922638E-01 
4.922638E-01 
4 . 922638E-01 
4 . 922638E-01 

2 

4.777190E-01 
4 . 833899E-01 
4 . 836864E-01 
4. 837045E-01 
4 . 837054E-01 
4 . 837056E-01 
4 . 837057E-01 
4 . 837058E-01 
4 . 837058E-01 

3 

4 . 576990E-01 
4 . 647904E-01 
4 . 651281E-01 
4 . 651474E-01 
4 . 651484E-01 
4 . 651487E-01 
4 . 651489E-01 
4.651490E-01 
4.651491E-01 

4 

3 . 972311E-01 
4 . 035690E-01 
4 . 039019E-01 
4 . 039205E-01 
4 . 039216E-01 
4 . 039218E-01 
4 . 039220E-01 
4 . 039220E-01 
4. 039220E-01 

5 

3.448500E-01 
3 . 501960E-01 
3 . 504831E-01 
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7 . 824449E-02 
7.824496E-02 
7 . 824505E-02 
7 . 824509E-02 
7 . 824510E-02 
7 . 824510E-02 

SLAB 

4. 083015E-02 
4. 166200E-02 
4 . 170692E-02 
4. 170934E-02 
4. 170950E-02 
4 . 170954E-02 
4. 170955E-02 
4. 170956E-02 
4. 170956E-02 

File  o”2: 


3 . 504989E-01 
3. 504998E-01 
3 . 505001E-01 
3 . 505002E-01 
3 . 505002E-01 
3 . 505002E-01 

6 

2 . 913136E-01 
2 . 955713E-01 
2 . 957992E-01 
2 . 958117E-01 
2 . 958124E-01 
2 . 958126E-01 
2.958127E-01 
2 . 958128E-01 
2 . 958128E-01 


FN  SOLUTION 
BDRY  FLUXES 

SLAB  NR=  1  X* [  O.OOOOOOOOOOOOOOOOE+OO, 

1 . OOOOOE+OO  2 . 10001E-01  1.80043E-02 

9. OOOOOE-Ol  2 . 23886E-01  1.99501E-02 

8.00000E-01  2 . 39473E-01  2.23673E-02 

7 . OOOOOE-Ol  2 . 56974E-0 1  2 . 54508E-02 

6 . OOOOOE-Ol  2.76546E-01  2.95198E-02 

5. OOOOOE-Ol  2 . 98152E-01  3.51361E-02 

4. OOOOOE-Ol  3 . 21268E-01  4.33882E-02 

3. OOOOOE-Ol  3.44242E-01  5.66942E-02 

2. OOOOOE-Ol  3 . 63308E-01  8.17101E-02 

1. OOOOOE-Ol  3.72672E-01  1.45459E-01 

l.OOOOOE-lO  3 . 59389E-01  3.74860E-01 

SLAB  NR=  2  X® [  5 . OOOOOOOOOOOOOOOOE-02 , 

1. OOOOOE+OO  2.01847E-01  3.55405E-02 

9. OOOOOE-Ol  2 . 15593E-01  3.92761E-02 

8. OOOOOE-Ol  2 . 31 1 17E-01  4.38882E-02 

7. OOOOOE-Ol  2 . 48678E-0 1  4.97254E-02 

6. OOOOOE-Ol  2 . 6851 IE-01  5.73492E-02 

5. OOOOOE-Ol  2 . 90703E-01  6.77247E-02 


5 . OOOOOOOOOOOOOOOOE-02] 


O.IOOOOOOOOOOOOOOO  ] 
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4 . 00000E-01  3.14918E-01  8.26587E-02 

3 . 00000E-01  3.39769E-01  1.05961E-01 

2 . 00000E-01  3.61746E-01  1.47170E-01 

1.00000E-01  3 . 75371E-01  2.36852E-01 

1.00000E-10  3.74861E-01  3.79362E-01 

SLAB  NR*  3  X* [  0.1000000000000000  ,  0.2000000000000000  ] 

1 . 00000E+00  1 . 92840E-01  6.82461E-02 

9 . 00000E-01  2. 06342E-01  7.50249E-02 

8.00000E-01  2 . 21674E-01  8.32905E-02 

7.00000E-01  2 . 39138E-01  9.35877E-02 

6 . 00000E-01  2.59037E-01  1.06761E-01 

5 . 00000E-01  2.81573E-01  1.24186E-01 

4 . 00000E-01  3.06585E-01  1.48252E-01 

3 . 00000E-01  3.32928E-01  1.83410E-01 

2 . 00000E-01  3 . 57283E-01  2.38446E-01 

1.00000E-01  3 . 74077E-01  3.26739E-01 

1 . 00000E-10  3.79340E-01  3.77342E-01 

SLAB  NR*  4  X* [  0.2000000000000000  ,  0.5000000000000000  ] 

1 . OOOOOE+OO  1.73234E-01  1.43159E-01 

9 . 00000E-01  1.86021E-01  1.55001E-01 

8 . 00000E-01  2 . 00692E-01  1.68882E-01 

7 . 00000E-01  2 . 17620E-01  1.85332E-01 

6 . OOOOOE-Ol  2 . 37232E-01  2.05036E-01 

5 . 00000E-01  2.59941E-01  2.28845E-01 

4. 00000E-01  2.85937E-01  2.57648E-01 

3. OOOOOE-Ol  3 . 14584E-01  2.91658E-01 

2. OOOOOE-Ol  3.42985E-01  3.27321E-01 

1. OOOOOE-Ol  3.64933E-01  3.47041E-01 

1.00000E-10  3 . 77347E-01  3.32752E-01 

SLAB  NR*  5  X* [  0.5000000000000000  ,  0.7500000000000000  ] 

1. OOOOOE+OO  1.08052E-01  1.78448E-01 

9. OOOOOE-Ol  1.17326E-01  1.90788E-01 

8. OOOOOE-Ol  1 . 28290E-01  2.04711E-01 

7. OOOOOE-Ol  1 . 41428E-01  2.20418E-01 

6. OOOOOE-Ol  1 . 57408E-01  2.38029E-01 

5. OOOOOE-Ol  1.77150E-01  2.57410E-01 

4. OOOOOE-Ol  2.01883E-01  2.77711E-01 

3. OOOOOE-Ol  2.32986E-01  2.96256E-01 

2. OOOOOE-Ol  2.70683E-01  3.06246E-01 

1. OOOOOE-Ol  3 . 08158E-01  2.96869E-01 

1.00000E-10  3.32752E-01  2.72345E-01 

SLAB  NR*  6  X*[  0 . 7500000000000000  ,  1.000000000000000  ] 
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1 . OOOOOE+OO 
9.00000E-01 
8.00000E-01 
7 . 00000E-01 
6 . 00000E-01 
5.00000E-01 
4.00000E-01 
3.00000E-01 
2 . OOOOOE-Ol 
1 . 00000E-01 
1.00000E-10 


5.20546E-02 
5.71218E-02 
6.32738E-02 
7.08964E-02 
8.05795E-02 
9.32668E-02 
1 . 10553E-01 
1.35283E-01 
1.72635E-01 
2.28599E-01 
2.72345E-01 


1.90265E-01 
2.00705E-01 
2 . 11884E-01 
2. 23634E-01 
2 . 35537E-01 
2.46679E-01 
2 . 55203E-01 
2.57582E-01 
2.48328E-01 
2 . 24033E-01 
1.86138E-01 


File  o3: 

SLAB  NR=  1  X=[  O.OOOOOOOOOOOOOOOOE+OO,  5 . OOOOOOOOOOOOOOOOE-02] 

0 . OOOOOOOOOOOOOOOOE+OO 
1. OOOOOE+OO  2 . 10001E-01  0. OOOOOE+OO 

9. OOOOOE-Ol  2 . 23886E-01  0. OOOOOE+OO 

8. OOOOOE-Ol  2 . 39473E-01  0. OOOOOE+OO 

7. OOOOOE-Ol  2.56974E-01  0. OOOOOE+OO 

6. OOOOOE-Ol  2 . 76546E-01  0. OOOOOE+OO 

5. OOOOOE-Ol  2.98152E-01  0. OOOOOE+OO 

4. OOOOOE-Ol  3 . 21268E-01  0. OOOOOE+OO 

3. OOOOOE-Ol  3.44242E-01  0. OOOOOE+OO 

2. OOOOOE-Ol  3.63308E-01  0. OOOOOE+OO 

1. OOOOOE-Ol  3.72672E-01  0. OOOCOE+OO 

1 . 00000E-10  3.59389E-01  0. OOOOOE+OO 

SLAB  NR=  2  X=[  5 . 0000000000000000E-02 ,  0.1000000000000000  ] 

5 . OOOOOOOOOOOOOOOOE-02 
1. OOOOOE+OO  2 . 01847E-01  1.80043E-02 

9. OOOOOE-Ol  2.15593E-01  1.99501E-02 

8. OOOOOE-Ol  2 . 31 1 17E-01  2.23674E-02 

7. OOOOOE-Ol  2.48678E-01  2.54508E-02 

6. OOOOOE-Ol  2.68511E-01  2.95198E-02 

5. OOOOOE-Ol  2.90703E-01  3.51361E-02 

4. OOOOOE-Ol  3 . 14918E-01  4.33882E-02 

3. OOOOOE-Ol  3 . 39769E-01  5.66943E-02 

2. OOOOOE-Ol  3.61746E-01  8.17103E-02 

1. OOOOOE-Ol  3 .7537 IE-01  1.45459E-01 

1 . OOOOOE-IO  3 . 74861E-01  3.74860E-01 

SLAB  NR=  3  X=[  0.1000000000000000  ,  0.2000000000000000  ] 
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0.1000000000000000 

1 . 00000E+00  1 . 92840E-01  3.55405E-02 

9 . OOOOOE-Ol  2 . 06342E-01  3.92762E-02 

8 . 00000E-01  2 . 21674E-01  4.38883E-02 

7. 00000E-01  2 . 39138E-01  4.97255E-02 

6. OOOOOE-Ol  2 . 59037E-01  5.73493E-02 

5. OOOOOE-Ol  2.81573E-01  6.77247E-02 

4. OOOOOE-Ol  3.06585E-01  8.26587E-02 

3. OOOOOE-Ol  3.32928E-01  1.05961E-01 

2. OOOOOE-Ol  3 . 57283E-01  1.47170E-01 

1. OOOOOE-Ol  3.74077E-01  2.36851E-01 

1.00000E-10  3.79340E-01  3.79361E-01 

SLAB  NR®  4  X=[  0.2000000000000000  ,  0.5000000000000000 

0.2000000000000000 

1 . OOOOOE+OO  1 . 73234E-01  6.82461E-02 

9. OOOOOE-Ol  1.86021E-01  7.50249E-02 

8. OOOOOE-Ol  2 . 00692E-01  8.32905E-02 

7. OOOOOE-Ol  2 . 17620E-01  9.35877E-02 

6. OOOOOE-Ol  2 . 37232E-01  1.06761E-01 

5. OOOOOE-Ol  2 . 59941E-01  1.24186E-01 

4. OOOOOE-Ol  2.85937E-01  1.48252E-01 

3. OOOOOE-Ol  3 . 14584E-01  1.83410E-01 

2. OOOOOE-Ol  3 . 42985E-01  2.38446E-01 

1. OOOOOE-Ol  3 . 64933E-01  3.26739E-01 

1.00000E-10  3 . 77347E-01  3.77341E-01 

SLAB  NR®  5  X® [  0 . 5000000000000000  ,  0.7500000000000000 

0.5000000000000000 

1. OOOOOE+OO  1 . 08052E-01  1.43159E-01 

9. OOOOOE-Ol  1.17326E-01  1.55001E-01 

8. OOOOOE-Ol  1 . 28290E-01  1.68882E-01 

7. OOOOOE-Ol  1.41428E-01  1.85332E-01 

6. OOOOOE-Ol  1 . 57408E-01  2.05036E-01 

5. OOOOOE-Ol  1 . 77 150E-01  2.28845E-01 

4. OOOOOE-Ol  2 . 01883E-01  2.57648E-01 

3. OOOOOE-Ol  2 . 32986E-01  2.91658E-01 

2. OOOOOE-Ol  2 . 70683E-01  3.27321E-01 

1. OOOOOE-Ol  3.08158E-01  3.47041E-01 

1 . 0000QE-10  3 . 32752E-01  3.32751E-01 

SLAB  NR*  6  X® [  0.7500000000000000  ,  1.000000000000000 

0.7500000000000000 

1. OOOOOE+OO  5 . 20546E-02  1.78448E-01 

9. OOOOOE-Ol  5.71218E-02  1.90788E-01 


8.00000E-01 

6.32738E-02 

2.04711E-01 

7.00000E-01 

7.08964E-02 

2.20418E-01 

6 . 00000E-01 

8 . 05795E-02 

2.38029E-01 

5 .00000E-01 

9 . 32668E-02 

2 . 57410E-01 

4 . 00000E-01 

1 . 10553E-01 

2.77711E-01 

3.00000E-01 

1.35283E-01 

2 . 96256E-01 

2 . 00000E-01 

1 . 72635E-01 

3 . 06246E-01 

1 .00000E-01 

2 . 28599E-01 

2 . 96869E-01 

1 .OOOOOE-IO 

2.72345E-01 

2 . 72345E-01 

1.000000000000000 

1 .00000E+00 

0 . 00000E+00 

1 . 90265E-01 

9 . 00000E-01 

0 . 00000E+00 

2 . 00705E-01 

8.00000E-01 

0 . 00000E+00 

2. 11884E-01 

7.00000E-01 

0 . 00000E+00 

2 . 23634E-01 

6 . 00000E-01 

0 . OOOOOE+OO 

2 . 35537E-01 

5 . 00000E-01 

0 . 00000E+00 

2 . 46679E-01 

4 . 00000E-01 

0. 00000E+00 

2 . 55203E-01 

3 . 00000E-01 

0.  OOOOOE+OO 

2 . 57582E-01 

2 . 00000E-01 

0. OOOOOE+OO 

2 . 48328E-01 

1 .00000E-01 

0.  OOOOOE+OO 

2 . 24033E-01 

1 .00000E-10 

0. OOOOOE+OO 

1 .86138E-01 
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Flux  versus  Direction  for  MultiSlab 
FN  Method 


-1.00  -0.73  -0.30  -0-23  0.00  0.25  0.50  0.75  1.00 

Direction  (mu) 


Figure  1:  Regeneration  of  Siewert’s  Data  from  Sample  Problem  Data 
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V.D  Detailed  Program  Notes 

i 

Detailed  programming  notes  are  given  in  this  section.  These  notes  refer  to  comments  and 
detailed  instructions  deemed  necessary  for  the  programmer  who  would  like  to  understand 
and  modify  the  MSlab  code. 


V.D.l  Main  Program 

These  are  the  comments  for  the  Main  Program. 


1.  CmainI  ***  Read  Input 

(a)  The  following  are  the  current  limits  contained  in  PARAMETER  statements  at  the 
beginning  of  the  main  program  and  where  required  in  subroutines: 

NP  =  Maximum  number  of  N  of  F;v  approximation  (70). 

ILM  =  Maximum  order  of  Gauss-Legendre  (G-L)  quadrature  (100). 

IMM1  =  Maximum  number  of  angular  edit  points  (100). 

INS  =  Maximum  number  of  slabs  (10). 

NPX  =  Maximum  number  of  spatial  edit  points  within  each  slab  (50). 

The  program  dimensions  are  automatically  adjusted  by  changing  the  values  of 
the  above  variables  in  the  PARAMETER  statements. 

(b)  If  only  one  slab  is  being  considered  (NS  =1),  the  number  of  boundary  iterations 
( LIT)  is  set  to  1. 

(c)  The  G-L  quadrature  order  required  to  determine  the  boundary  or  interior  fluxes 
can  be  different  for  each  slab  but  the  G-L  quadrature  order  for  the  evaluation  of 
integrals  over  the  boundary  fluxes  connecting  the  slabs  will  always  be  the  same 
for  each  slab. 

(d)  Four  possibilities  for  basis  functions  have  been  provided.  Table  S  lists  these 
functions.  From  numerical  experiment,  the  value  of  7  (XAR)  =  0.75  (for  IBR 
=  3  or  4)  has  been  found  to  give  the  best  results. 

2.  CmainS  ***  F.v  Iteration  Loop 

(a)  This  is  the  main  iteration  loop  which  cycles  through  the  F.v  approximation  foi  j 

increasing  N. 

(b)  N  is  incremented  by  L12D  specified  on  Line  2  of  the  input. 
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3.  CmainT  ***  Boundary  Flux  Iteration  Loop 

Since  the  boundary  fluxes  at  the  slab  interfaces  for  the  multiple  slab  case  are  not 
known  a  priori ,  an  iterative  process  is  required  for  their  determination.  By  numerical 
experiment,  it  has  been  verified  that  rather  than  having  the  F/v  iteration  serve  also  as 
a  boundary  flux  iteration,  it  is  more  efficient  to  include  at  least  several  boundary  flux 
iterations  at  each  N.  This  reduces  the  N  required  for  convergence.  However,  there 
is  a  tradeoff  since  if  many  boundary  flux  iterations  are  included,  the  computational 
time  will  be  greatly  increased.  In  addition  the  advantage  this  iteration  provides  is 
reduced  as  N  increases. 

4.  CmainS  ***  Slab  Loop 

This  loop  cycles  through  the  slabs  with  the  previously  determined  boundary  fluxes. 

5.  CmainIO  ***  Check  To  Avoid  Same  Quadrature  Abscissa  Edit  Points  and 
Collocation  Points 

Since  0.5  is  a  zero  of  all  odd  order  shifted  Legendre  polynomials,  a  modification 
must  be  made  to  avoid  division  by  zero  in  some  parts  of  the  calculation.  This  is 
accomplished  by  not  allowing  the  F/v  order  N  (NN)  and  the  quadrature  order  for  the 
determination  of  the  matrix  elements  (LMR)  to  both  be  odd  (or  even)  simultaneously. 
If  they  are  both  odd  or  even,  then  LMR  is  increased  by  one.  A  similar  procedure  is 
required  for  the  quadrature  order  for  the  integration  of  the  boundary  fluxes  (LM3). 

C.  Cmain13  ***  Specify  End  or  Interior  Slab  and  Boundary  Conditions  (bc) 

In  subroutine  Si,  a  slab  interface  indicator  for  the  right  (IBDR)  and  left  (IBDL) 
interfaces  is  set  according  to  the  values  in  Table  9.  These  flags  are  necessary  for  the 
appropriate  integration  of  the  boundary  fluxes  (subroutines  INTO,  INTR,  INTRZ) 
to  form  the  inhomogeneous  terms  of  the  collocation  equations. 

7.  CmainIG  ***  Find  Matrix  Elements  and  Inhomogeneous  Terms 

When  IFLG  =  1,  the  matrix  elements  and  inhomogeneous  terms  are  formed  in  the 
symmetric  form  in  subroutine  MELB  to  be  used  in  the  boundary  flux  determination 
[eqs.  (25)  and  eqs.  ( 26 ) j -  For  IFLG  =  -1,  the  matrix  elements  and  inhomogeneous 
terms  are  put  in  the  asymmetric  form  for  eqs.  (20). 

8.  CmainIO  ***  Check  F,v  Convergence  For  Boundary  Fluxes  at  the  Edit 
Points 

Before  the  calculation  of  the  interior  flux  is  init  iated,  the  boundary  fluxes  must  have 
converged  relative  to  the  last  F/v  determination  at  the  desired  angular  edit  point. 
Once  the  boundary  fluxes  have  converged,  then  the  interior  fluxes  are  determined 
until  they  converge  at  all  angular  and  spatial  edit  points.  Even  though  the  boundary 
fluxes  have  converged,  they  are  recalculated  at  the  same  F/v  approximation  for  the 
interior  fluxes  since  they  are  required  at  the  same  collocation  points. 
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Table  8:  Basis  Functions  Available  to  Program  and  How  to  Access  Them 


IB 

Name 

'I’c(n) 

n 

Shifted  Legendre  polynomials 

Pa(2/i  -  1) 

2 

Monomials 

3 

Shifted  Monomials 

WBsam 

4 

Modified  Legendre  polynomials 

Pa(2/^-l) 

Table  9:  Interface/Boundary  Condition  and  How  to  Access  Them 


IBDR  or  IBDL 

Interface/Boundary  Condition 

0 

Slab /Slab  no  source 

1 

Vacuum/Slab  with  beam  source 

2 

Slab/ Vacuum  no  source 

V.D.2  FNX  Subprogram 


These  are  the  comments  for  the  FNX  Subprogram.  Its  main  purpose  is  to  calculate  the 
values  of  the  flux  in  the  interior  of  a  slab. 


1.  CfnxI  ***  Initialize  Convergence  Flag  Ips 

CfnxI  ***  Ips(Ir+Nr,Ix)  =  -1  Point  Ix  in  Slab  Nr  Not  Yet  Converged 
CfnxI  ***  Ips(Ir+Nr,Ix)  =  1  Point  Ix  in  Slab  Nr  has  Converged 

Initialize  all  convergence  pointers  to  the  not  converged  value  (-1). 

2.  Cfnx2  ***  Determine  Matrix  Elements  for  Current  Collocation  Points 
and  Required  Sums 

When  IFLG  =  -2,  the  matrix  elements  are  calculated  in  subprogram  MELB.  The  A 
and  B  function  values  are  found  and  summed. 

3.  Cfnx3  ***  Spatial  Grid  Loop 

If  the  point  has  already  converged,  then  it  is  skipped.  If  not,  then  determine  where 
it  is  in  the  slab  and  set  up  some  flags  for  the  next  step. 

4.  Cfnx4  ***  Set  Up  Matrix  Elements  and  Inhomogeneous  Term  for  Inte¬ 
rior  Flux  Calculation 

For  IFLG  =  -1,  the  matrix  elements  and  inhomogeneous  terms  are  put  in  the  asym¬ 
metric  form  of  eqs.  (20)  and  calculated  using  the  MELB  subprogram. 

5.  Cfnx5  ***  Solve  for  the  Coupling  Coefficients 

Performs  the  matrix  inversion  required  to  solve  for  the  coupling  coefficients  in  eqs.  (20). 

6.  Cfnx6  ***  Determine  Interior  Fluxes 

Using  the  post-processing  algorithm,  the  interior  fluxes  are  found  using  eqs.  (21). 

7.  CfnxT  ***  Test  for  Convergence 

If  the  relative  error  from  the  last  iteration  is  below  the  tolerance,  then  set  the  con¬ 
verged  flag.  If  not.  then  continue. 
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ABSTRACT 

A  set  of  benchmark  calculations  is  described  for  the  one-dimensional,  inultigioup 
discrete  ordinates  method  as  applied  to  electron  transport.  An  analytical  representation 
of  the  group  scalat  flux  is  obtained  by  application  of  the  Fourier  transform  operator  to 
the  transport  equation,  expansion  of  the  transformed  equation  in  terms  of  the  flux 
moments,  and  re- inversion  of  the  transformed  scalar  flux  moments.  The  results  of  S„ 
calculations,  with  both  the  lineal  discontinuous  and  diamond  difference  algorithms,  are 
compared  with  benchmark  calculations  for  a  plane  isotropic  source  of  200  keV  electrons 
in  infinite  aluminum.  The  degtee  of  scattering  anisotropy  is  allowed  to  vary  from 
isotropic  to  third-order  anisotropic.  Analysis  of  spatial  discretization  error  is  also 
discussed. 

I.  Introduction 

Ail  integral  part  of  transport  method  development  is  veiification  that  the  numerical 
algorithm  used,  as  well  as  its  coding,  perform  piopcrly.  Unfortunately,  due  to  the  lack 
of  appropriate  bcnchmaiks  and  the  fact  that  verification  efforts  are  usually  tangential  to 
the  main  thrust  of  the  development,  code  verification  is  often  not  seriously  considered. 
Of  course,  a  full  veiification  is  not  possible  since  if  it  were,  there  would  be  no  need  for 
the  method  development  in  the  fir st  place.  Partial  verification  however,  can  be 
achieved  by  application  of  any  oi  all  of  the  following  procedures: 

a)  monitoring  of  particle  conservation; 

b)  comparison  with  intuitive  ticnds; 

c)  comparison  to  similar  algorithms; 

d)  comparison  to  analytical  solutions. 

It  should  be  emphasized  that  the  success  of  a  verification  procedure  does  not  provide 
proof  positive  that  the  method  is  operational.  However,  the  lack  thereof  does  indicate 
that  the  numerical  algorithm  or  its  coding  is  >n  some  way  defective  Thi‘  piesentation 
is  concerned  with  the  establishment  of  a  benchmark  for  1-D  multigroup  tiansport 
algorithms  and  the  verification  of  a  specific  S„  algorithm  u«ed  in  electron  transport 
theory,  as  an  example 

II.  Mulligroup  Benchmark  Formulation 

Tin'  slowing  down  and  deflection  of  elections  in  a  medium  larger  than  the  election 
range  (effectively  a  planar  infinite  medium)  can  be  approximated  hv  'lie  following 
inultigioup  formulation,  where  all  the  symbols  have  their  usual  meaning:1  1 
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The  group  mean  free  path  Xg  and  group  transfer  cross  sections  E((.  arc  determined 
assuming  that  continuous  slowing  down  theory  is  applicable1.  A  source  has  been  placed 
at  the  center  of  the  infinite  medium  and  is  assumed  to  emit  electrons  in  the  highest 
energy  group  with  a  beam  or  isotropic  angular  distribution 
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The  analytical  representation  of  the  group  scalar  flux  *•  (x)  is  obtained  by 
application  of  the  Fourier  transform  operator  to  Eq.  (1)  in  the  form 


and  its  subsequent  inversion 


,(k.„)  =  |dNe-ikxye(x.,)  , 

-oc 

no 

,(x.„)  =  £jdk  eikxv,(k,p). 


In  (he  transform  k-space,  Eq.(l)  becomes 
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where 
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Eq  (3)  can  be  iccast  in  the  following  form  for  the  moments  : 
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The  solution  to  Eqs  (3)  in  matrix  form  is 

«,=  A,-’q9, 
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where  the  elements  of  6g ,  Ag  and  q9  are  respectively, 

#Jil-(2f+I)StMT»,(k). 

and  ’£  ^(2C+1)E(  .  TJ.,(k)  3,.,((k)  +  Q,«#  ,  . 
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The  desired  scalar  flux  is  therefore 


Of  course,  one  could  attempt  to  evaluate  Eq  (C)  by  complex  contour  methods,  but 
for  a  cross  section  dependent  on  energy,  this  would  be  a  hopeless  task  for  all  but  a  few 
group  calculation.  For  this  reason,  a  direct  numerical  implementation  of  the  inversion 
is  performed.  Briefly,  evaluation  of  Eq.(6)  makes  use  of  the  equivalent  form 
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The  integral  is  then  expressed  as  an  infinite  series  with  each  term  being  an  integration 
over  cither  the  half  or  full  periods  of  the  cosine.  Each  integral  is  numerically  performed 
using  an  efficient  Romberg  integrator.  Finally,  the  series  is  numerically  summed  via 
the  Euler- Kuopp  transformation  to  accelerate  convergence.3 

III.  S„  Method 

The  S„  code  used  in  this  study  was  written  to  serve  primarily  as  a  general  purpose 
slab-geometry  research  code  for  electron  transport  studies,  and  as  a  test-bed  for 
multigroup  algorithms.  It  incorporates  the  same  linear  discontinuous  and  diamond 
difference  algorithms  as  are  found  in  ONETRAN4.  For  simplicity,  the  S„  code  was 
written  to  apply  only  to  slab  geometry.  The  main  departure  of  our  code  from 
ONETRAN  is  that  the  angular  components,  rather  than  the  Legendre  moments,  of  the 
flux  and  source  and  angular  scattering  kernels  are  retained  throughout  the  entire 
calculation.  This  feature  can  enable  straightforward  incorporation  of  angularly 
dependent  total  cross  sections  and  source  terms  which  may  be  used  to  represent  the 
presence  of,  for  example,  an  externally  applied  electric  field5.  This  code  is  currently 
being  used  in  two  areas  of  election  transport  research:  1)  calculation  of  ionization  doses 
resulting  from  the  incidence  of  gamma  radiation  on  multilayered  material  media  as  may 
be  typified  by  microelectronic  devices;  and  2)  studies  of  low  energy  electron-phonon 
scattering  in  solids. 

IV.  Results 

The  quantity  taken  as  the  basis  for  comparison  of  the  S„  and  licnrhmark 
calculations  is  the  scalar  flux  The  appioach  taken  in  the  comparisons  between 
niultigioup  S„  and  benchmark  calculations  begins  with  a  sciies  of  one- group  calculations 
foi  a  plane  isotropic  souicc  imbedded  in  an  essentially  infinite  aluminum  medium  in 
which  the  scattering  cross  section  was  taken  to  be  linearly  anisotropic.  For  this  one- 
group  case  we  investigated  the  behavior  of  the  ratio  of  successive  errors  as  the  spatial 
mesh  size  was  decreased  in  older  to  reproduce  the  known  theory.  Then  for  the 
multigroup  cases.  S„  and  benchmark  calculations  were  made  for  a  plane  isotropic  source 
imbedded  in  an  infinite  aluminum  medium.  These  calculations,  carried  out  using  40 
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uniformly  spaced  electron  energy  groups,  were  made  for  a  200  keV  electron  source  in 
aluminum.  Four  scattering  cases,  corresponding  to  progressively  higher  scattering 
anisotropy  order,  were  investigated.  In  the  first,  an  artificial  isotropic  scattering  cross 
section  was  devised,  and  in  the  remaining  four,  S„  and  benchmark  calculations  were 
compared  for  extended  transport  corrected  cross  sections  of  order  2  through  4. 

Because  the  benchmark  calculations  can  only  accommodate  an  infinite  niedium, 
and  the  S„  algorithm  requires  a  finite  niedium,  we  chose  the  dimensions  of  the  test 

mesh  portions  of  the  scattering  media  to  he  deeply  imbedded  within  two  very  thick 

coarsely  meshed  slabs  each  of  thickness  10  mean  free  paths.  For  the  multigroup 
calculations,  we  chose  the  slab  dimensions  to  correspond  to  ±1  r.u.  (  =  range  unit), 
which  for  our  200  keV  clcctrou  source  in  aluminum  translates  into  a  thickness  of 
±  0.058  g/cmJ.  These  geometries  closely  approximate  infinite  media  for  all  practical 
purposes. 

One-Group  Calculation: 

Plane  Isotropic  Source;  Linearly  Anisotropic  Scattering 

A  set  of  5  S,j  calculations  were  made  with  the  linear  discontinuous  differencing 
algorithm  for  a  plane  isotropic  somcc  located  at  the  center  (x=0)  of  a  thick  scattering 
medium  extending  20  mean  free  paths  in  both  positive  and  negative  x-dircctions.  The 

scattering  parameters  were  (  in  units  of  mean  free  path),  L'(  =  10,  ~  (18,  t,  = 

Zs  1 3.  The  spatial  mesh  discietization  step  sizes  for  the  5  calculations  wer?  taken  tci  W 
iix„  =  2J-";  11=1.2.3,4.5.  In  each  calculation.  Ax  was  held  constant  over  the  spatial 
mesh  test  region  (-10  <x  <  410).  Figure  1  is  a  plot  of  the  successive  error  ratios  r„(x) 
(  =  »„(x)/»  .  ,(x)J.  From  the  curves  shown,  a  region  of  uniformity(3  <  |x|  <  G)  could  be 
indcntified  lor  successive  discietization  error  ratios.  In  this  region  it  was  found  that  the 
error  always  decreased  for  n  =  1.2.3.  For  the  n  =  l  and  n=2  cases,  the  ratio  is  nearly 
constant  and  has  the  value  ~  8,  which  indicates  third  order  accurate  behavior  for  the 
linear  discontinuous  algorithm.  Also  in  this  region,  when  n=3.  it  is  seen  that  the  error 
ratio  settles  down  to  —  S  for  one-half  of  the  interval;  however,  the  ratio  r4  does  not  show 
a  consistent  decrease  in  error  as  Ax  is  decreased  from  0.25  to  0.125  This  is  attributable 
to  the  fact  that  it  is  not  possible  to  further  improve  the  result  by  halving  Ax  with  a 
fixed  inner  iteration  convergence  precision.  The  fractional  errors  observed  in  these  cases 
were  generally  comparable  to  the  convergence  precision  ( 10  -  s).  Large  variations  in  i„ 
were  observed  both  near  the  source  plane  and  the  region  boundary.  The  poor  results 
near  the  source  plane  arise  from  the  singular  nature  of  the  source  geometry,  while  near 
the  region  boundary,  changes  in  mesh  size  propagate  errors  into  the  test  zone.  The 
uniform  region  is  situated  well  away  from  both. 

Multigroup  Calculations: 

Plane  Isotropic  Source;  Isotropic  Scattering 

Four  sets  of  S„  calculations  were  made  foi  the  200  keV  isotropic  election  source 
imbedded  in  infinite  aluminum.  The  scnttciing  cross  section  employed  for  these 
calculations  is  physically  unrealistic  in  that  200  keV  electrons  in  aluminum  do  not 
scatter  isotropically.  However  in  this  first  attempt  to  test  the  multicroup  benchmark, 
we  chose  for  convenience  a  cross  vetion  which  contains  only  the  zero-th  order  Legendre 
coefficient  of  the  P-12  extended  transport  conectrd  Mott  cross  section  (subsequent 
calculations,  repotted  in  the  next  section,  make  use  of  the  actual  P-2.3.4  extended 
transport  corrected  cross  sections).  These  calculations  covered  the  following  four  cases: 
linear  discontinuous  and  diamond  differencing,  each  used  with  S6  and  Su  Gauss- 
Legendre  quadrature  sets.  In  each  of  these  four  cases.  5  separate  S„  calculations  were 
made.  corres;>onding  to  a  graduated  scale  of  spatial  mesh  resolutions  |ax  =  02.  .01, 
.004.  .002,  .001  r.u  ]  spanning  a  distance  of  0  2  r.u.  fiom  the  source  plane  at  x=0.  The 
reasons  for  the  restriction  to  0.2  r.u.  are  twofold:  1)  to  eliminate  vacuum  boundary 


effects  in  the  S„  calculations;  2)  to  operate  within  the  high  accuracy  range  of  the 
analytical  benchmark  solution. 

In  the  following  we  show  only  S6  results,  since  in  all  cases  studied  the  difference 
between  S6  and  Su  was  found  to  be  negligible.  Figure  2  is  a  plot  of  the  scalar  flux  as 
calculated  by  the  benchmark  method  (solid  lines)  and  the  S„  rode  (circles)  with  the 
diamond  difference  scheme  for  spatial  resolution  Ax=.004  r.u.  The  10  curves  are  the 
energy  spectra  at  x  =  .02(.U2).2  r.u.  As  can  be  seen,  the  agreement  is  excellent. 

Figure  3  is  a  plot  of  the  exponent,  b,  of  the  relative  error  dependence  |(  a.  (Ax)"1] 
of  the  spatial  resolution  vs.  group  midpoint  energy,  for  6  positions  along  x,  as  the 
resolution  is  increased  fiom  Ax  =  .02  to  Ax=  0!  r.u.  for  the  linear  discontinuous 
algorithm.  The  relative  error  (<)  is  evaluated  as  the  deviation  of  the  S„  result  from  the 
benchmark  result.  The  inner  iteration  convergence  precision  for  the  upper  graph  was 
set  to  10  "3,  while  for  the  lower  graph,  we  chose  10"  .  As  anticipated,  the  exponent  is 
-3  as  dictated  by  theory.  It  snould  be  noted  that  this  holds  true  for  a  wider  energy 
range  as  the  accuracy  of  the  S„  calculation  is  increased.  However,  at  some  point 
sufficiently  low  on  the  energy  scale,  no  further  improvement  is  achieved.  As  might  be 
expected,  the  more  accurate  S„  calculation,  to  the  extent  that  this  crossover  behavior 
exists  at  all,  exhibits  this  behavior  at  a  lower  energy  value.  We  believe  that  this  is 
explained  by  an  accumulation  of  the  downscattcr  source  error,  which  is  not  predominant 
in  the  high  energy  groups.  It  can  also  be  noted  that  tightening  of  the  convergence 
precision  mitigates  most  of  the  cnatic  behavior  of  the  exponent  at  the  two  positions 
(x  =  .02,  .06  r.u.)  closest  to  the  souice  plane. 

Plane  Isotropic  Source;  Anisotropic  Scattering 

The  200  keV  electron  souice  (at  x=0)  and  aluminum  scatterer  configuration 
(|x|  <  1.0  r.u.)  chosen  for  the  anisotropic  scattering  benchmark  tests  were  the  same  as 
tliose  used  in  the  previous  case.  It  has  already  been  demonstrated6- 1  that  a  realistic 
description  of  the  transport,  or  at  least  the  resultant  energy  deposition,  of  200  keV 
electrons  in  aluminum  can  be  achieved  using  a  P-12  extended  transport  corrected 
screened-Ruthcrford  scattering  cioss  section  in  conjunction  with  the  S„  method. 
Doubtless,  the  use  of  lower  order  transport  corrections  for  this  same  problem  will  results 
in  a  sacrifice  of  physical  realism.  However  successive  application  of  the  extended 
transport  correction  from  P-1  to  P-4  solves  the  ‘'same"  problem  if  only  very 
approximately,  and  additionally  provides  us  with  a  means  for  systematically  increasing 
the  scattering  anisotropy  order  of  the  benchmark  calculations.  Figures  4-6  show  plots 
of  the  scalar  flux  at  10  equispaced  x  values  [0.04(0.04)0.4]  for  the  P-2,  P-3  and  P-4 
extended  transport  corrected  cross  sections  (the  L-th  Legendre  coefficient  of  the  P-L 
extruded  transport  cioss  section  is  zero,  hence  these  represent  scattering  anisotropy 
oideis  1,2.3).  In  all  cases  the  linear  discontinuous  differencing  algorithm  was  employed 
in  40  group-Su  calculations.  The  spatial  discretization  step  was  held  constant  at 
Ax=0.002  r.u.  As  before,  the  solid  curves  represent  the  benchmark  solution.  The 
agreement  between  the  two  methods  ranges  between  the  third  and  fifth  significant 
figure  for  all  points  except  at  the  souice  plane. 

V.  Concluding  Remarks 

We  believe  that  this  investigation  has  served  as  a  dcmonstiation  of  the  first 
multigioup  bcnchmatk  transport  solution.  The  excellent  agreement  between 
multigroup  S„  and  the  benchmark  lesults  provides  verification  of  the  predicted  error 
behavior  for  both  the  linear  discontinuous  and  diamond  difference  schemes. 
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Fig.  1.  Ratios  of  successive  errors  for  5  one-group  linear  discontinuous  S,j 
calculations  Spatial  discretization  step  size  is  varied  according  to  Ax„=2:~n\ 
n  =  l-5.  Stable  region,  3  <  |x|  <  6  (mean  free  path  units),  is  indicated. 
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Fig.  3  Plots  of  the  exponent  (61  of  the  relative  error  [»],  where  (  »  ix"‘,  as  -ix  is 
decreased  from  .02  to  .01  r.u.  The  behavior  of  b  is  shown  as  a  function  of  energy 
at  the  six  x  positions  indicated  for  two  linear  discontinuous  Sg  calculations  with 
inner  iteration  convergence  precisions  given  by  10_,(upper  graph)  and  10"5(lower 
graph). 
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Fig  4.  Comparison  of  the  scalar  flux,  F*(E)  vs.  E,  for  40  energy  groups  at  10  x 
positions  as  indicated  Solid  lines  represent  the  benchmark  calculation  results; 
circles  represent  the  results  of  a  linear  discontiuous  S,j  calculation.  Anisotropic 
scattering  of  200  keV  electrons  in  A1  with  P-2  extended  transport  corrected  cross 
sections. 
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Fig.  5.  Comparison  of  the  scaiar  flux,  FA(E)  vs.  E,  for  40  energy  groups  at  10  x 
positions  as  indicated.  Solid  lines  represent  the  benchmark  calculation  results; 
circles  represent  the  results  of  a  linear  discontiuous  calculation  Anisotropic 
scattering  of  200  keV  electrons  in  A1  with  P-3  extended  transport  corrected  cross 
sections. 
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Fig.  6.  Comparison  of  the  scalar  flux,  F*(E)  vs.  E,  for  40  energy  groups  at  10  x 
positions  as  indicated.  Solid  lines  represent  the  benchmark  calculation  results; 
circles  represent  the  results  of  a  linear  discontiuous  Sjj  calculation.  Anisotropic 
scattering  of  200  keY  electrons  in  A!  with  P-4  extended  transport  corrected  cross 
sections. 
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